海洋捕食者算法 MPA

前言

这个算法介绍很少，所以想边学习边记录一下。可能更新有点慢（近期大概率不更新），各位小伙伴们见谅。
原论文: Marine Predators Algorithm: A nature-inspired metaheuristic
代码: https://github.com/afshinfaramarzi/Marine-Predators-Algorithm

本文链接: 原创
https://blog.csdn.net/weixin_43850253/article/details/110724958

思想

海洋捕食者算法启发于自然界中捕食者的捕食策略。该算法认为顶级捕食者具有最大的搜索本领。
顶级捕食者们构成精英矩阵(一个顶级捕食者即为问题的一个解)。

海洋捕食者算法（MPA）是一种新型元启发式优化算法。参考作者论文, 给出算法的流程:

(1) 初始化精英矩阵(Elite)和猎物矩阵(Prey)

猎物矩阵(Prey) 矩阵每一个元素 X_ij 的初始化方法:
$$X_{ij}=X_{min}+rand(X_{max}-x_{min})$$
最终得到的Prey矩阵:

其中，n是种群的规模，d是每个维度的位置(问题的解的维度)。

对每一个Prey个体X_i = [X_i,1, X_i,2, …, X_i,d], 计算其适应度, 然后使用适应度最优的个体 X^I 复制n份构成Elite矩阵

其中n是种群的规模，d是每个维度的位置(问题的解的维度)，Elite的维度与Prey的维度相同。

(2)接着我们开始进行优化。在优化的过程中，具有三个步骤。

步骤一:

当迭代次数小于最大迭代次数的三分之一的时候,
$$\begin{gathered} s_i=R_B⨂(Elit{e}_i-R_B⨂Pre{y}_i),i=1...n \\ Pre{y}_i=Pre{y}_i+P.R⨂s_i \end{gathered}$$
其中，R_B 是采用布朗随机游走产生的随机数组成的向量，维度是 d(问题的求解规模，下同)。s_i 代表移动的步长。 P是一个常数，等于0.5。R是一个0到1之间的均匀分布的随机数组成的向量，维度是 d。

R_B相当于一般化的高斯分布(Normal Gaussian distribution)。每一个元素 R_Bi 可以通过下列表达式来计算:
$$R_{Bi}=\frac{1}{\sqrt{2\pi }}exp(-\frac{x^2}{2})$$

步骤二:

当迭代次数大于最大迭代次数的三分之一而小于其三分之二时，种群分两部分进行操作。

前半部分种群跟新规则如下：
$$\begin{gathered} s_i=R_L⨂(Elit{e}_i-R_L⨂Pre{y}_i),i=1,...,n/2 \\ Pre{y}_i=Pre{y}_i+P.R⨂s_i \end{gathered}$$
其中， R_L 是 Levy 分布组成的出来的一个向量，维度是 d。P是一个常数，等于0.5。R是一个0到1之间的均匀分布的随机数组成的向量，维度是 d。

R_L 的每一项元素 R_Li 可以由下列式子计算得来:
$$R_{Li}=C\times \frac{x}{y^{1/a}}$$
其中，C 和 α是一个常数，分别等于0.05和1.5。
$$\begin{gathered} x=Normal(0,{\sigma }_x^2) \\ y=Normal(0,{\sigma }_y^2) \end{gathered}$$
在上面的表达式中
$$\begin{gathered} {\sigma }_x=[\frac{\Gamma (1+\alpha )sin(\frac{\pi \alpha }{2})}{\Gamma (\frac{1+\alpha }{2})\alpha 2^{\frac{\alpha -1}{2}}}{]}^{1/\alpha } \\ {\sigma }_y=1 \\ \alpha =1.5 \end{gathered}$$

后半部分种群跟新规则如下:
$$\begin{gathered} s_i=R_B⨂(R_B⨂Elit{e}_i-Pre{y}_i),i=n/2,...,n \\ Pre{y}_i=Elit{e}_i+P.CF⨂s_i \end{gathered}$$
这里R_B 是采用布朗随机游走产生的随机数组成的向量，维度是 d。P是常数，等于0.5。CF是步长s_i 的自适应参数(下同), 定义为
$$CF=(1-\frac{Iter}{Max\_Iter}{)}^{(2\frac{Iter}{Max\_Iter})}$$
其中, Iter是迭代次数，Max_Iter是最大迭代次数。

步骤三：

当迭代次数大于最大迭代次数的三分之二时，进入第三个阶段，此时种群更新规则如下:
$$\begin{gathered} s_i=R_L⨂(R_L⨂Elit{e}_i-Pre{y}_i),i=1,...,n \\ Pre{y}_i=Elit{e}_i+P.CF⨂s_i \end{gathered}$$

(3) 解决涡流形成和FADS效应(Eddy formation and FADs’ effect)

此操作的作用是让算法在迭代过程中尽可能跳出局部最优解，已达到更好的寻优精度。
$$Pre{y}_i=\{\begin{array}{l}Pre{y}_i+CF[X_{min}+R⨂(X_{max}-X_{min})]⨂U,ifr<=FADS\\ Pre{y}_i+[FADs(1-r)+r](Pre{y}_{r1}-Pre{y}_{r2}),ifr>FADS\end{array}$$
其中r是一个随机数, FADS是一个影响优化过程的常数，等于0.2。r1和r2是Prey两个随机下标, 1 ≤ r1,r2 ≤ n。 U是一个包含0和1的二进制向量，维度是d。U的每一个元素 U_i 定义为
$$U_i=\{\begin{array}{l}0,ifrandom\le FADs\\ 1,ifrandom>FADs\end{array}$$
其中random是一个0到1的随机数，FADs等于0.2。

(4) 海洋记忆(Marine memory)

这一步骤进行对Elite(精英)矩阵的更新。

针对每一个Prey矩阵中的个体Prey_i ,计算其适应度，若适应度由于Elite矩阵矩阵中相应的位置的适应度时，则将该个体替代原来精英矩阵中相应的个体。

然后在计算整个精英矩阵中最优个体的适应度，若符合要求，则算法结束，否则继续迭代。

算法的流程图总结出来如下:

当然，Afshin Faramarzi[6]等人在论文中也提到，这只是MPA算法的第一个版本，后续还可以继续改进。尽管如此MPA在他们和其他学者的实验中都发现，MPA算法具有很高的寻优精确度。

李代华和崔东文[7]将海洋捕食者算法(MPA)与自适应神经模糊推理系统(ANFIS)相结和的方法来预测径流。Mohammed A. A. Al-qaness[8] 等人用MPA算法来预测意大利，美国，伊朗和韩国的新冠肺炎确诊病例。

代码

2022.4.21更，这是我2020年写的代码，当时为了应付课程作业，还没来得及优化和排小bug

# %%
import numpy as np
import pandas as pd  # 读取csv文件的库
import matplotlib.pyplot as plt
import random
import os

cur = os.getcwd().replace("\\", '/')
if cur.split('/')[-1] != 'MPA':
    os.chdir('D:/project/pyCharmProject/cal_intelligent/report/MPA')
print(os.getcwd())

# %%
filename = 'data/ENB_nor.xlsx'
if filename.split('.')[-1] == 'csv':
    df = pd.read_csv('../' + filename, encoding='gbk')
elif filename.split('.')[-1] == 'xlsx' or filename.split('.')[-1] == 'xls':
    df = pd.read_excel('../' + filename, sheet_name=0)
else:
    df = 'error! 不支持的文件格式'
print(df)

# %%
x_m = df[['X1', 'X2', 'X3', 'X4', 'X5', 'X6', 'X7', 'X8']].values.tolist()
y_lt = df['Y1'].tolist()
print('shape of x_m: {} * {}'.format(len(x_m), len(x_m[0])))
print('shape of y_lt: {}'.format(len(y_lt)))

# %%
length = len(x_m)
x = x_m[:int(length * 0.8)]
y = y_lt[:int(length * 0.8)]


# %%
def fc_MSE(w, x, y):
    """
    适应度函数
    :param w: w是一系列参数的组合, 其中, 0到31是ANN输入层到隐藏层的32个w权重，32到35是输入层到隐藏层的4个偏置b,
    36到39是隐藏层到输出层的4个w权重, 40是隐藏层到输出层的1个偏置b, 因此len(w) = 41
    :param x: 特征矩阵, n * 8 维度, n是数据条数
    :param y: 标签, n * 1 维度, n是数据条数
    :return: 适应度值
    """
    # 注意切片运算时左闭右开 w[a:b] -> [a,b)
    w1 = np.array(w[:32]).reshape((8, 4))
    b1 = np.array(w[32:36]).reshape((1, 4))
    w2 = np.array(w[36:40]).reshape((4, 1))
    b2 = np.array(w[40])
    x_n = np.array(x).reshape((len(x), 8))
    y_n = np.array(y).reshape((len(x), 1))
    temp = x_n.dot(w1) + b1
    temp = 1 / (1 + np.exp(-temp))
    result = y_n - temp.dot(w2) + b2
    MSE = np.mean(np.square(y - result))

    # print('MSE: ', MSE)
    # https://blog.csdn.net/qq_42257962/article/details/108265730
    # RMSE = np.sqrt(np.mean(np.square(y - y_hat)))
    # MAE = np.mean(np.abs(y - y_hat))
    # MAPE = np.mean(np.abs((y - y_hat) / y)) * 100
    return MSE


def Levy_lt(dimension: int):
    # gamma(1.25) = 0.9064
    # gamma(2.5) = 1.329340388179137
    sigma_x = 0.6965757716463266
    sigma_y = 1
    levy_list = []
    while True:
        levy_x = np.random.normal(loc=0, scale=sigma_x ** 2, size=1)[0]
        levy_y = np.random.normal(loc=0, scale=sigma_y ** 2, size=1)[0]

        levy_a = 0.05 * levy_x / (levy_y ** (1 / 1.5))
        if np.isnan(levy_a):
            continue
        levy_list.append(levy_a)
        if len(levy_list) >= dimension:
            break

    # print(levy_a)
    return np.array(levy_list)


# %%
# 这里用MPA来求解神经网络的参数, 求解的函数为使得 y - f(w, x) 的向量的MSE最小
def MPA(x, y, func: str):
    """
    使用 pso 来算出神经网络连接的初始化权重
    :param x: 特征矩阵, n * 8 维度, n是数据条数
    :param y: 标签, n * 1 维度, n是数据条数
    :param func: 适应度函数的函数名
    :return: 长度为41的 w, 将作为神经网络的权重
    """
    # 先定义一些参数
    x_min, x_max = -5, 5  # 解的范围, 及将要传入神经元的参数的范围
    N = 4  # 种群规模
    D = 41  # 觅食空间，解的维度
    Max_Iter = 300  # 最大迭代次数，建议500~1000
    FADs = 0.2  # the probability of FADs effect on the optimization process
    P = 0.5  # 在 Max_Iter/3 代数内中会用到

    # 开始操作
    # Prey 代表 Prey矩阵
    Prey = np.array([(x_min + np.random.random((D, 1)) * (x_max - x_min)).reshape(41).tolist()
                     for i in range(N)])
    # print(Prey)
    # print(len(Prey))  # 4
    # print(len(Prey[0]))  # 41

    # 存放精英矩阵中最好的解
    g_best = [i for i in Prey[0]]
    g_best_value = eval(func)(Prey[0], x, y)

    # 找出现阶段全局最优解
    for i in range(1, N):
        temp = eval(func)(Prey[i], x, y)
        if temp < g_best_value:
            g_best_value = temp
            g_best = [j for j in Prey[i]]

    # 再将最好的x_i复制 N 份构成精英矩阵
    Elite = np.array([[i for i in g_best] for j in range(N)])
    # for i in Elite:
    #     print(i)
    # 存放精英矩阵中每一行的解的值
    e_value_lt = [g_best_value for i in range(N)]

    print('g_best: {}\ng_best_value: {}'.format(g_best, g_best_value))
    print()

    generation_num = 0
    while generation_num < Max_Iter:
        generation_num += 1
        CF = (1 - generation_num / Max_Iter) ** (2 * generation_num / Max_Iter)

        # 开始三阶段演变
        # 第一阶段
        if generation_num < Max_Iter / 3:
            R_B = np.random.normal(loc=0, scale=1, size=D)
            R = np.random.random(size=D)
            for i in range(N):
                step_size = R_B * (Elite[i] - R_B * Prey[i])  # numpy中矩阵间*代表相同位置元素相乘
                Prey[i] += P * R * step_size

        # 第二阶段
        if Max_Iter / 3 <= generation_num < 2 * Max_Iter / 3:
            for i in range(int(N / 2)):
                R_L = Levy_lt(D)
                R = np.random.random(size=D)
                step_size = R_L * (Elite[i] - R_L * Prey[i])
                Prey[i] += P * R * step_size
            for i in range(int(N / 2), N):
                R_B = np.random.normal(loc=0, scale=1, size=D)
                step_size = R_B * (R_B * Elite[i] - Prey[i])
                Prey[i] = Elite[i] + P * CF * step_size

        # 第三阶段
        if generation_num >= 2 * Max_Iter / 3:
            for i in range(N):
                R_L = Levy_lt(D)
                step_size = R_L * (R_L * Elite[i] - Prey[i])
                Prey[i] = Elite[i] + P * CF * step_size

        # Eddy formation and FADs’ effect 增加跳出局部最优的能力
        random_num = random.random()
        if random_num <= FADs:
            U = np.array([0 if random.random() <= 0.2 else 1 for i in range(D)])
            for i in range(N):
                R = np.random.random(size=D)
                Prey[i] += CF * (x_min + R * (x_max - x_min)) * U
        else:  # random_num > FADs
            for i in range(N):
                temp = np.random.randint(0, N, size=2)
                # r1 and r2 subscripts denote random indexes of prey matrix.
                r1, r2 = temp[0], temp[1]
                Prey[i] += (FADs * (1 - random_num) + random_num) * (Prey[r1] - Prey[r2])

        # 评估和替换
        for i in range(N):
            temp = eval(func)(Prey[i], x, y)
            # 比表上相应位置还要优
            if temp < e_value_lt[i]:
                Elite[i] = [value for value in Prey[i]]
                e_value_lt[i] = temp
                # 比全局最优还要优
                if temp < g_best_value:
                    g_best = [value for value in Prey[i]]
                    g_best_value = temp
                    print('gen_num_{}更新啦...\n g_best: {}, g_best_value: {}'.format(generation_num, g_best,
                                                                                   g_best_value))
    print('final g_best: {}\ng_best_value: {}'.format(g_best, g_best_value))
    return g_best, g_best_value


MPA(x, y, "fc_MSE")

结语

参考文献: a paper:
Marine Predators Algorithm: A nature-inspired metaheuristic
Afshin Faramarzi…