1 算法思想

黏菌算法由李世民等人发表于2020年，模拟了黏菌觅食过程中的行为和形态变化。

黏菌在有丝分裂后形成的变形体成熟之后，进入营养生长时期，会形成网状型态，且依照食物、水与氧气等所需养分改变其表面积。在黏菌算法中，黏菌会根据当前位置的客观条件（适应度函数优劣），决定每个个体所在位置的权重，然后个体会根据权重决定新的位置在哪。

当黏菌接近食物源时，生物振荡器会通过静脉产生传播波，来增加细胞质流量。食物浓度越高，生物振荡器产生的传播波越强，细胞质流动越快。黏菌算法就是通过模拟黏菌这种捕食行为来实现智能寻优功能的。

借鉴黏菌的生物行为，可以抽象出三个规则：

1. 接近食物： 黏菌通过空气中的气味接近食物，黏菌接近食物时呈圆形与扇形结构运动。
2. 包围食物： 黏菌静脉接触的食物浓度越高，生物振荡器产生的传播波越强，细胞质流动越快。
3. 抓取食物： 黏菌在食物浓度低时更慢地接近食物，找到优质食物时更快地接近食物。

2 算法步骤

1. 设置参数，初始化种群，计算适应度值；
2. 更新黏菌权重W、参数a、参数b：
   
   W
   =
   {
   1
   +
   r
   ×
   l
   o
   g
   (
   b
   F
   −
   S
   (
   i
   )
   b
   F
   −
   w
   F
   +
   1
   )
   ,
   i
   f
   该
   个
   体
   的
   适
   应
   度
   值
   排
   在
   群
   体
   前
   一
   半
   1
   −
   r
   ×
   l
   o
   g
   (
   b
   F
   −
   S
   (
   i
   )
   b
   F
   −
   w
   F
   +
   1
   )
   ,
   e
   l
   s
   e
   在
   后
   一
   半
   W=\left\{ \begin{aligned} 1+r×log(\frac{bF-S(i)}{bF-wF}+1) , & if 该个体的适应度值排在群体前一半 \\ 1-r×log(\frac{bF-S(i)}{bF-wF}+1) , & else 在后一半 \end{aligned} \right.
   W=⎩⎪⎪⎨⎪⎪⎧​1+r×log(bF−wFbF−S(i)​+1),1−r×log(bF−wFbF−S(i)​+1),​if该个体的适应度值排在群体前一半else在后一半​
   （下面有参数说明表）
   公式解读：前一半较优部分个体的权重取值较大，在[1,1.3]内；较差的后一半取值在[0.7,1]内，越差的个体权重越接近0.7。取log的作用就是减缓了W的变化率。
   
   a
   =
   a
   r
   c
   t
   a
   n
   h
   (
   b
   )
   b
   =
   1
   −
   t
   T
   a=arctanh(b) \\ b=1-\frac{t}{T}
   a=arctanh(b)b=1−Tt​
   b的变化是一条线性递减的简单直线，a的变化如下图：
   
3. 更新个体位置$X(t+1)$，分三种情况：
   
   X
   (
   t
   +
   1
   )
   =
   {
   r
   ×
   (
   u
   b
   −
   l
   b
   )
   +
   l
   b
   ,
   r
   <
   z
   X
   b
   (
   t
   )
   +
   v
   b
   (
   W
   ×
   X
   A
   (
   t
   )
   −
   X
   B
   (
   t
   )
   )
   ,
   z
   ≤
   r
   <
   p
   v
   c
   ×
   X
   (
   t
   )
   ,
   p
   ≤
   r
   X(t+1)=\left\{ \begin{aligned} &r×(ub-lb)+lb, & r<z \\ &X_b(t)+vb(W×X_A(t)-X_B(t)), & z \le r<p \\ &vc×X(t), & p \le r \end{aligned} \right.
   X(t+1)=⎩⎪⎨⎪⎧​​r×(ub−lb)+lb,Xb​(t)+vb(W×XA​(t)−XB​(t)),vc×X(t),​r<zz≤r<pp≤r​
   （下面有参数说明表）
   若r<z，按第一个子公式更新位置；否则更新$p,vb,vc$，并比较r和p，p的计算如下：
   $$p=tanh(|S(i)-DF|),i=1,2,...,Nvb=[-a,a]，vc=[-b,b]$$
   p的变化如下图：
   
   若r<p，按第二个子公式更新位置；否则，按第三个子公式更新位置。

更新个体位置这一步模拟了黏菌的生物行为：利用生物振荡器产生的传播波改变静脉中的细胞质流动速度。
算法是这样模拟的：通过vb，vc，W来模拟静脉宽度的变化和振荡器振荡频率变化，当食物浓度低时，慢慢接近食物，扩大全局搜索能力，当找到优质食物时，迅速接近食物，加强局部搜索能力。
公式解读：第一个子公式获取的是全局随机位置，类似于GA中的变异操作；第二个是在当前最优位置的附近搜索，类似于局部搜索；第三个没看懂，好像会让个体最优值能收敛到0，当最优解不为0时效果不好。

4. 计算适应度值，更新全局最优解；
5. 判断是否满足结束条件，若不满足，返回步骤2。

参数说明表：

3 求函数最值（Python实现）

求解下列函数的最小值：

f
(
x
1
,
x
2
)
=
x
1
2
+
x
2
2
f(x_1,x_2)=x_1^2+x_2^2
f(x1​,x2​)=x12​+x22​
标准答案是0，函数长这个样子：

主函数：

import numpy as np
from matplotlib import pyplot as plt
import SMA
'''适应度函数'''
def fun(X):
    Results = np.sum(X ** 2)
    return Results
'''主函数 '''
# 设置参数
pop = 30  # 种群数量
MaxIter = 500  # 最大迭代次数
dim = 2  # 维度
lb = -10 * np.ones(dim)  # 下边界
ub = 10 * np.ones(dim)  # 上边界
# 调用SMA算法
GbestScore, GbestPositon, Curve = SMA.SMA(pop, dim, lb, ub, MaxIter, fun)
print('最优适应度值：', GbestScore)
print('最优解[x1,x2]：', GbestPositon)
# 绘制适应度曲线
plt.figure(1)
plt.plot(Curve, 'r-', linewidth=2)
plt.xlabel('Iteration', fontsize='medium')
plt.ylabel("Fitness", fontsize='medium')
plt.grid()
plt.title('SMA', fontsize='large')
plt.show()

SMA.py：

import numpy as np
import copy as copy
'''黏菌优化算法'''
'''
Args:
    pop: 种群数量
    dim: 个体维度
    lb: 下边界，维度[1,dim]
    ub: 上边界，维度[1,dim]
    MaxIter: 最大迭代次数
    fun: 适应度函数接口
Returns:
    GbestScore: 最优解对应的适应度值
    GbestPositon: 最优解
    Curve: 画迭代曲线用的
'''
def SMA(pop, dim, lb, ub, MaxIter, fun):
    # 1.设置参数，初始化种群，并计算适应度值
    # 设置参数
    z = 0.03  # 一个参数，更新个体位置时用，表示随机分布的黏菌个体占总体的比例
    Curve = np.zeros([MaxIter, 1]) # 画迭代图用
    W = np.zeros([pop, dim])  # 黏菌权重
    # 初始化种群
    X = initialization(pop, ub, lb, dim)
    # 计算适应度值
    fitness = CaculateFitness(X, fun)  # 计算适应度值
    fitness, sortIndex = SortFitness(fitness)  # 对适应度值排序
    X = SortPosition(X, sortIndex)  # 种群排序
    GbestScore = copy.copy(fitness[0])
    GbestPositon = copy.copy(X[0, :])
    for t in range(MaxIter):
        worstFitness = fitness[-1]
        bestFitness = fitness[0]
        S = bestFitness - worstFitness + 10E-8  # 当前最优与最差适应度的差值，算权重W用，10E-8为极小值，避免分母为0；
        # 2.更新权重W和参数a、参数b
        # 权重W
        for i in range(pop):
            if i < pop / 2:  # 适应度值排前一半的W计算
                W[i, :] = 1 + np.random.random([1, dim]) * np.log10((bestFitness - fitness[i]) / (S) + 1 + 10E-8)
            else:  # 适应度值排后一半的W计算
                W[i, :] = 1 - np.random.random([1, dim]) * np.log10((bestFitness - fitness[i]) / (S) + 1 + 10E-8)
        # 参数b
        b = 1 - (t / MaxIter)
        # 参数a
        if b != -1 and b != 1:
            a = np.math.atanh(b)
        else:
            a = 1
        # 3.位置更新，分三种情况
        for i in range(pop):
            # 3.1如果r<z，按第一个子公式更新位置
            if np.random.random() < z:
                X[i, :] = (ub.T - lb.T) * np.random.random([1, dim]) + lb.T
            # 否则，更新参数p，vb，vc，继续判断
            else:
                p = np.tanh(abs(fitness[i] - GbestScore))
                vb = 2 * a * np.random.random([1, dim]) - a
                vc = 2 * b * np.random.random([1, dim]) - b
                for j in range(dim):
                    r = np.random.random()
                    A = np.random.randint(pop) # 随机选择两个黏菌个体
                    B = np.random.randint(pop)
                    # 3.2如果r<p，按第二个子公式更新位置
                    if r < p:
                        X[i, j] = GbestPositon[j] + vb[0, j] * (W[i, j] * X[A, j] - X[B, j])
                    # 3.3按第三个子公式更新位置
                    else:
                        X[i, j] = vc[0, j] * X[i, j]
        X = BorderCheck(X, ub, lb, pop, dim)  # 边界检查
        fitness = CaculateFitness(X, fun)  # 计算适应度值
        fitness, sortIndex = SortFitness(fitness)  # 对适应度值排序，得到排好序的适应度值和对应的索引
        X = SortPosition(X, sortIndex)  # 根据排好序的索引对种群排序
        if (fitness[0] <= GbestScore):  # 更新全局最优
            GbestScore = copy.copy(fitness[0])
            GbestPositon = copy.copy(X[0, :])
        Curve[t] = GbestScore
    return GbestScore, GbestPositon, Curve

运行结果：
最优适应度值： [6.28823104e-226]
最优解： [[-1.77578646e-113 1.77054045e-113]]
可以看到答案非常接近最优适应度值0。

4 算法进阶

集思广益：

1. W、vc、vb等参数的设置可以使个体在任意角度形成搜索向量，即在任意方向上搜索解空间，使算法具有找到最优解的可能性。
2. 文中对于帮助算法及时跳出局部最优值并没有做出很多的工作。
3. CEC2017效果一般。
4. CEC2014前几个函数效果很好，比别的算法可以优于几十个数量级，观察了一下这几个函数的最优值都是0，猜测可能和那个位置更新公式刻意逼近0有关。
5. CEC2014剩余函数效果一般，说明探索能力和跳出局部最优能力有待改善。
6. 黏菌算法没有贪心步骤，无论新位置好不好，个体都移动到新位置。

直接改进SMA

融合别的智能优化算法来改进SMA

SMA及其改进的应用

参考书籍：范旭，《Python智能优化算法——从原理到代码实现与应用》第一版，电子工业出版社。