黏菌算法（Slime Mould Algorithm，SMA）

1 算法思想

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

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

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

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

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

2 算法步骤

1. 设置参数，初始化种群，计算适应度值；
2. 更新黏菌权重W、参数a、参数b：
   $$W=\{\begin{array}{rl}1+r\times log(\frac{bF-S(i)}{bF-wF}+1),& if该个体的适应度值排在群体前一半\\ 1-r\times log(\frac{bF-S(i)}{bF-wF}+1),& else在后一半\end{array}$$
   （下面有参数说明表）
   公式解读：前一半较优部分个体的权重取值较大，在[1,1.3]内；较差的后一半取值在[0.7,1]内，越差的个体权重越接近0.7。取log的作用就是减缓了W的变化率。
   $$\begin{gathered} a=arctanh(b) \\ b=1-\frac{t}{T} \end{gathered}$$
   b的变化是一条线性递减的简单直线，a的变化如下图：
   
3. 更新个体位置$X(t+1)$，分三种情况：
   $$X(t+1)=\{\begin{array}{rlr}& r\times (ub-lb)+lb,& r<z\\ & X_b(t)+vb(W\times X_A(t)-X_B(t)),& z\le r<p\\ & vc\times X(t),& p\le r\end{array}$$
   （下面有参数说明表）
   若r<z，按第一个子公式更新位置；否则更新$p,vb,vc$，并比较r和p，p的计算如下：
   $$\begin{gathered} p=tanh(|S(i)-DF|),i=1,2,...,N \\ vb=[-a,a]，vc=[-b,b] \end{gathered}$$
   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$$
标准答案是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()

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
29
30
31

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

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
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85

运行结果：
最优适应度值： [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智能优化算法——从原理到代码实现与应用》第一版，电子工业出版社。