人工智能|海洋捕食者算法 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) 矩阵每一个元素 Xij 的初始化方法:
X i j = X m i n+r a n d ( X m a x?x m i n ) X_{ij} = X_{min} ~+~ rand(X_{max} ~-~ x_{min}) Xij?=Xmin? + rand(Xmax? ? xmin?)
最终得到的Prey矩阵:

文章图片
其中，n是种群的规模，d是每个维度的位置(问题的解的维度)。
对每一个Prey个体Xi = [Xi,1, Xi,2, …, Xi,d], 计算其适应度, 然后使用适应度最优的个体 XI 复制n份构成Elite矩阵

文章图片
其中n是种群的规模，d是每个维度的位置(问题的解的维度)，Elite的维度与Prey的维度相同。
(2)接着我们开始进行优化。在优化的过程中，具有三个步骤。
步骤一:
当迭代次数小于最大迭代次数的三分之一的时候,
s i = R B ?( E l i t e i ? R B ? P r e y i ) ,i = 1... n P r e y i = P r e y i+P . R ? s is_i = R_B \bigotimes ~(Elite_i - R_B \bigotimes Prey_i), ~i = 1...n \\ Prey_i = Prey_i ~+~ P.R \bigotimes s_i ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ si?=RB?? (Elitei??RB??Preyi?), i=1...nPreyi?=Preyi? + P.R?si?
其中，RB 是采用布朗随机游走产生的随机数组成的向量，维度是 d(问题的求解规模，下同)。si 代表移动的步长。 P是一个常数，等于0.5。R是一个0到1之间的均匀分布的随机数组成的向量，维度是 d。
RB相当于一般化的高斯分布(Normal Gaussian distribution)。每一个元素 RBi 可以通过下列表达式来计算:
R B i = 1 2 π e x p ( ? x 2 2 ) R_{Bi} = \frac{1}{\sqrt{2\pi}}exp(-\frac{x^2}{2}) RBi?=2π ?1?exp(?2x2?)
步骤二:
当迭代次数大于最大迭代次数的三分之一而小于其三分之二时，种群分两部分进行操作。
前半部分种群跟新规则如下：
s i = R L ?( E l i t e i?R L ? P r e y i ) ,i = 1 , . . . , n / 2 P r e y i = P r e y i+P . R ? s is_i = R_L \bigotimes ~(Elite_i ~-~ R_L \bigotimes Prey_i),~i=1,...,n/2 \\ Prey_i = Prey_i ~+~ P.R \bigotimes s_i ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ si?=RL?? (Elitei? ? RL??Preyi?), i=1,...,n/2Preyi?=Preyi? + P.R?si?
其中， RL 是 Levy 分布组成的出来的一个向量，维度是 d。P是一个常数，等于0.5。R是一个0到1之间的均匀分布的随机数组成的向量，维度是 d。
RL 的每一项元素 RLi 可以由下列式子计算得来:
R L i = C×x y 1 / a R_{Li} = C ~\times~ \frac{x}{y^{1/a}} RLi?=C × y1/ax?
其中，C 和 α是一个常数，分别等于0.05和1.5。
x = N o r m a l (0 ,σ x 2) y = N o r m a l (0 ,σ y 2) x = Normal(~0,~ \sigma_x^2~) \\ y = Normal(~0,~ \sigma_y^2~) x=Normal( 0, σx2? )y=Normal( 0, σy2? )
在上面的表达式中
σ x = [ Γ ( 1 + α ) s i n ( π α 2 ) Γ ( 1 + α 2 ) α 2 α ? 1 2 ] 1 / α σ y = 1α = 1.5\sigma_x = [\frac{\Gamma(1 + α)sin(\frac{\piα}{2})}{\Gamma(\frac{1+α}{2})α2^{\frac{α-1}{2}}}]^{1/α} \\ \sigma_y = 1 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\\ α = 1.5 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ σx?=[Γ(21+α?)α22α?1?Γ(1+α)sin(2πα?)?]1/ασy?=1α=1.5
后半部分种群跟新规则如下:
s i = R B ?( R B ? E l i t e i ? P r e y i ) ,i = n / 2 , . . . , n P r e y i = E l i t e i + P . C F ? s is_i = R_B \bigotimes~(R_B \bigotimes Elite_i - Prey_i),~i=n/2,...,n \\ Prey_i = Elite_i + P.CF \bigotimes s_i ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ si?=RB?? (RB??Elitei??Preyi?), i=n/2,...,nPreyi?=Elitei?+P.CF?si?
这里RB 是采用布朗随机游走产生的随机数组成的向量，维度是 d。P是常数，等于0.5。CF是步长si 的自适应参数(下同), 定义为
C F = ( 1 ? I t e r M a x _ I t e r ) ( 2 I t e r M a x _ I t e r ) CF = (1 - \frac{Iter}{Max\_Iter})^{(2\frac{Iter}{Max\_Iter})} CF=(1?Max_IterIter?)(2Max_IterIter?)
其中, Iter是迭代次数，Max_Iter是最大迭代次数。
步骤三：
当迭代次数大于最大迭代次数的三分之二时，进入第三个阶段，此时种群更新规则如下:
s i = R L ?( R L ? E l i t e i?P r e y i ) , i = 1 , . . . , n P r e y i = E l i t e i + P . C F ? s is_i = R_L \bigotimes ~(R_L \bigotimes Elite_i ~-~ Prey_i), i=1,...,n \\ Prey_i = Elite_i + P.CF \bigotimes s_i~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ si?=RL?? (RL??Elitei? ? Preyi?),i=1,...,nPreyi?=Elitei?+P.CF?si?
(3) 解决涡流形成和FADS效应(Eddy formation and FADs’ effect)
此操作的作用是让算法在迭代过程中尽可能跳出局部最优解，已达到更好的寻优精度。
P r e y i = { P r e y i + C F [ X m i n + R ? ( X m a x ? X m i n ) ] ? U , i fr < = F A D S P r e y i + [ F A D s ( 1 ? r ) + r ] ( P r e y r 1 ? P r e y r 2 ) ,i f r> F A D S Prey_i = \left\{ \begin{array}{l} Prey_i + CF[X_{min} + R \bigotimes(X_{max} - X_{min})] \bigotimes U, if~ r <= FADS \\ Prey_i + [FADs(1-r)+r](Prey_{r1}-Prey_{r2}), ~if r~> FADS \end{array} \right. Preyi?={Preyi?+CF[Xmin?+R?(Xmax??Xmin?)]?U,if r<=FADSPreyi?+[FADs(1?r)+r](Preyr1??Preyr2?), ifr >FADS?
其中r是一个随机数, FADS是一个影响优化过程的常数，等于0.2。r1和r2是Prey两个随机下标, 1 ≤ r1,r2 ≤ n。 U是一个包含0和1的二进制向量，维度是d。U的每一个元素 Ui 定义为
U i = { 0 ,i fr a n d o m ≤ F A D s 1 ,i fr a n d o m > F A D s U_i = \left\{ \begin{array}{l} 0, ~~~if~ random ≤ FADs \\ 1, ~~~if~ random > FADs \end{array} \right. Ui?={0,if random≤FADs1,if random>FADs?
其中random是一个0到1的随机数，FADs等于0.2。
(4) 海洋记忆(Marine memory)
这一步骤进行对Elite(精英)矩阵的更新。
针对每一个Prey矩阵中的个体Preyi ,计算其适应度，若适应度由于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 oscur = 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 MSEdef 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 = https://www.it610.com/article/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 = https://www.it610.com/article/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_valueMPA(x, y, "fc_MSE")

结语