2、算法原理分析：

        HO算法使用了一个三阶段模型，结合了河马在河流或池塘中的位置更新、对捕食者的防御策略和逃避方法。该算法通过自适应地调整搜索空间的分辨率和搜索速度，以快速而准确地找到最优解，具有收敛速度快、求解精度高等特点。

        2.1 随机初始解生成

        在HO算法中，河马是优化问题的候选解，这意味着每个河马在搜索空间中的位置更新表示决策变量的值。因此，每一头河马都被表示为一个向量，而河马种群在数学上由一个矩阵来表征。类似于传统的优化算法，HO的初始化阶段涉及随机初始解的生成。在此步骤中，使用以下公式生成决策变量的向量：

 

其中Xi表示第i头河马在每个决策变量下的所有位置信息，为一个行向量。xij表示第i头河马在第j个决策变量下的位置信息，xij即表示一个候选解。r是0到1范围内的随机数，并且Ibj和ubj分别表示第j个决策变量的下限和上限。N表示种群中的河马数量，m表示问题中决策变量的数量（即x、y、z分别表示一个决策变量，求函数h=x^2+y^2+z^2的最小值时m=3）。

此时的初始解xij均为 [ lbj, ubj ] 范围内的随机数。（例如，在求y=x^2在[-10,10]内的最小值时。m=1（只有一个变量x）, lb = -10, ub = 10，河马群内河马的总头数为N。）

河马种群位置信息的矩阵如下：

        2.2 第一阶段（探索阶段）：河马在河川或池塘中的位置更新

        河马群由几只成年雌性河马、幼河马、多只成年雄性河马和占主导地位的雄性河马（河马群的领导者）组成。基于目标函数值迭代（最小化问题的最低值和最大化问题的最高值）来确定优势种群。一般来说，河马倾向于彼此靠近聚集。占优势的雄性河马保护河马群和领地免受潜在的威胁。多只雌性河马被安置在雄性河马周围。在达到成熟期后，雄性河马会被占统治地位的雄性河马从群体中驱逐出去。随后，这些被驱逐的雄性个体被要求要么吸引雌性，要么与河马群中其他已建立的雄性成员进行统治地位的竞争，以建立自己的统治地位。

以下方程表示了湖泊或池塘中雄性河马的位置信息：

代表雄性河马的位置，y1是0和1之间的随机数，代表具有支配地位河马的位置（在当前迭代中具有最佳函数值的河马），其中的等于1或2。

此时雄性河马的位置即为当前河马向具有支配地位河马的位置不同程度靠近的结果。（该过程可能会陷入局部最优解）

（若初始解中优势解为2，而当前解为4，则计算得到的当前雄性河马的位置为 4 + y1×(2 - I1×4)，即4再向2靠近。）

        以下方程表示了湖泊或池塘中雌性或未成熟河马的位置信息：

其中的h1和h2是从h方程中的五种情况中随机选择的数或向量。I2等于1或2。MGi指的是一些从河马群中随机选择的河马的平均值，这些河马被选择的概率相等，MGi则表示了。

r1到r4是0和1之间的随机向量，r5是0和1之间的随机数。h场景使用了I1和I2，增强了算法的全局搜索能力，提高了算法的探索性。它提供了更好的全局搜索能力。r7是0和1之间的随机数。Q1和Q2是整数随机数，可以是1或0。t表示当前迭代次数，而表示最高迭代次数。

        上式描述了雌性或未成熟河马在河马群中的位置（）。大多数未成熟的河马都在它们的母亲附近，但由于好奇心，有时未成熟的河马会与河马群分开远离它们的母亲。如果T大于0.6，这意味着未成熟的河马已经远离了它的母亲。r6是介于0和1之间的数，如果r6大于0.5，则表明未成年河马已远离其母亲，但仍在河马群内或附近，否则，它已脱离河马群。

（即当刚开始迭代时，T比较大，雌性或未成熟河马在解空间中进行了较大的探索，远离了当前的最优解或其局部邻域。随着迭代次数的增加T不断减小，当 r6​>0.5时，表明未成年河马虽然远离了母亲，但仍然在河马群内或附近，这在算法中对应于在解空间中进行中等程度的探索，即河马在保持与当前最优解一定距离的同时，也在探索新的区域。如果r6​≤0.5，则表示河马在解空间的边界内随机移动，模拟了河马完全脱离群体，进行更大范围的探索。）

        下式描述了雄性和雌性或未成熟河马在河马群中的位置更新。Fi为目标函数值。当雄性和雌性或未成熟河马的适应度值小于当前的适应度值时，当前的适应度值便进行更新。

        2.3 第二阶段（探索阶段）：河马防御捕食者

        河马群居生活的一个关键原因可以归结为它们的安全保障。这些庞大而沉重的动物群的存在可以阻止捕食者接近它们。然而，由于它们天生的好奇心，未成年的河马偶尔会偏离族群，成为尼罗河鳄鱼、狮子和斑点鬣狗的潜在目标，因为它们的力量相对于成年河马来说相对较弱。生病的河马和未成熟的河马一样，也容易被捕食者捕食。河马的主要防御策略是迅速转向捕食者，发出响亮的叫声，以阻止捕食者靠近它们。在这一阶段，河马可能会表现出接近捕食者的行为，以诱导其撤退，从而有效地抵御潜在的威胁。

下式表示捕食者在搜索空间中的位置：

其中表示从0到1的随机向量。

表示第i条河马到捕食者的距离。在这段时间里，河马会采取基于因子 的防御行为来保护自己免受捕食者的攻击。如果 小于，表明捕食者非常接近河马，在这种情况下，河马迅速转向捕食者并向其移动以使其撤退。如果 更大，则表明捕食者离河马的领地距离远，在这种情况下，河马转向捕食者，但运动范围更有限。其目的是让捕食者或入侵者意识到其在其领土内的存在。

        是河马面对捕食者时的位置。f是在2和4之间的均匀随机数，c是在1和1.5之间的均匀随机数，d是在2和3之间的均匀随机数。g表示在-1和1之间的均匀随机数。是维数为1 ×m的随机向量。是一个具有Levy分布的随机向量，用于在攻击河马时捕食者位置的突然变化。Levy运动的随机运动的数学模型计算如下：

w和v是位于[0,1]的随机数，是一个常数（ = 1.5），注意区分开和。是Gamma函数的缩写，可以通过下式得到：

下式为面对捕食者的河马在河马群中的位置更新：

如果，则意味着这只河马已经被猎杀，另一只河马将在兽群中取代它，否则捕食者将会逃跑，这只河马将回到兽群中。在第二阶段，全局搜索进程得到显著加强。第一阶段和第二阶段相互补充，有效地降低了陷入局部最小值的风险。

（捕食者的初始位置和河马一样是在解空间内的任意值。1、如果当前捕食者的适应度值比当前河马的适应度值小，则捕食者对河马的威胁较大，在这种情况下，河马会采取一种积极的应对策略，使用D的倒数来调整向捕食者移动行为的强度，这导致河马在解空间中进行较大的跳跃，以快速远离当前的局部最优解；2、如果当前捕食者的适应度值比当前河马的适应度值大，则捕食者对河马的威胁较小，在这种情况下，河马会采取一种更为保守的应对策略，通过增加一个随机扰动项 r9并结合 2×D 来调整向捕食者移动行为的强度，这导致河马在解空间中进行中等程度的探索，以保持与当前最优解一定距离的同时，也在探索新的区域。）

        2.3 第三阶段（开发阶段）：河马逃离捕食者

        河马面对捕食者的另一种行为是当河马遇到一群捕食者或者无法用它的防御行为击退捕食者时，河马会试图离开这个区域。通常情况下，河马会试图跑到最近的湖泊或池塘，以避免捕食者的伤害，因为狮子和鬣狗会避免进入湖泊或池塘。这种策略会使河马找到一个接近其当前位置的安全位置。在HO算法的第三阶段中建模这种行为，增强了算法在本地的搜索能力。为了模拟这种行为，在河马的当前位置附近生成随机位置。

是河马逃离捕食者时的位置，它被搜索以找到最近的安全位置。s1是从三种s情况中随机选择的随机向量或数。s所考虑的场景促使了更好的局部搜索。代表0到1之间的随机向量，而r10和r13表示在0和1的范围内产生的随机数。另外，r12是正态分布的随机数。t表示当前迭代次数，而表示最高迭代次数。

当新的位置提高了适应度值时，表明河马在其当前位置附近找到了更安全的位置，并相应地改变了其位置。

        在更新种群的HO算法中，没有将种群分成未成熟、雌性和雄性河马这三个单独的类别，因为虽然将它们分成单独的类别可以更好地模拟它们的性质，但这会降低优化算法的性能。在完成HO算法的每次迭代之后，基于阶段一至三更新所有群体成员，一直持续到最后一次迭代。在算法执行期间，最佳潜在解被一直地跟踪和存储，在整个算法完成后，最佳候选者即为问题的最终解决方案。

        探索（Exploration）指的是在解空间中寻找新的可能解，而开发（Exploitation）则是利用已知的最优解来进一步优化。河马优化算法通过模拟河马的行为，试图在探索和开发之间找到平衡。

        HO的过程细节的流程图和伪代码在下图中给出：

3、应用与代码：

论文地址：

• Hippopotamus optimization algorithm: a novel nature-inspired optimization algorithm | Scientific Reports

HO算法的matlab源代码可在MathWorks上公开获取，论文原作者已开源。

通过网盘分享的文件：HO_Matlab.zip
链接: https://pan.baidu.com/s/1MdF-kjnxhaKM9urnM-6HqA?pwd=feuv 提取码: feuv

HO算法的python代码，本人自己根据matlab代码和对算法的理解写了一个，如下：

import numpy as np
import random
import matplotlib.pyplot as plt
from scipy.special import gamma


####################################### 定义目标函数 #####################################################
def fun_info(F):
    if F == 'F1':
        fitness = F1
        lowerbound = -10
        upperbound = 10
        dimension = 30
    elif F == 'F2':
        fitness = F2
        lowerbound = -10
        upperbound = 10
        dimension = 30
    elif F == 'F3':
        fitness = F3
        lowerbound = -4
        upperbound = 4
        dimension = 1
    elif F == 'F4':
        fitness = F4
        lowerbound = -100
        upperbound = 100
        dimension = 2
    else:
        raise ValueError("Unknown function F specified")

    return lowerbound, upperbound, dimension, fitness


# F1 function
def F1(x):
    return np.sum(x ** 2)


# F2 function
def F2(x):
    return np.sum(np.abs(x)) + np.prod(np.abs(x))

def F3(x):
    return np.sin(x) + np.cos(x)

def F4(X):
    x = X[0]
    y = X[1]
    result = x**2 + y**2 + 10
    return result

###################################### end #####################################################


####################################### HO函数 #####################################################
def levy(n, m, beta):
    # 计算分子
    num = gamma(1 + beta) * np.sin(np.pi * beta / 2)
    # 计算分母
    den = gamma((1 + beta) / 2) * beta * 2 ** ((beta - 1) / 2)
    # 计算标准偏差
    sigma_u = (num / den) ** (1 / beta)
    # 生成正态分布的随机变量 u
    u = np.random.normal(0, sigma_u, (n, m))
    # 生成正态分布的随机变量 v
    v = np.random.normal(0, 1, (n, m))
    # 计算 levy 步长 z
    z = u / (np.abs(v) ** (1 / beta))
    return z


def HO(search_agents, max_iterations, lowerbound, upperbound, dimension, fitness_func):

    # 初始化 X 矩阵，大小为 (SearchAgents, dimension)
    X = np.zeros((search_agents, dimension))
    X_P1 = np.zeros((search_agents, dimension))
    X_P2 = np.zeros((search_agents, dimension))
    X_P3 = np.zeros((search_agents, dimension))
    X_P4 = np.zeros((search_agents, dimension))

    # 生成 X 矩阵中的值
    for i in range(dimension):
        # 对于每个搜索代理，生成一个维度为 dimension 的随机数组
        # 每个元素都在 [lowerbound, upperbound] 范围内
        X[:, i] = lowerbound + np.random.rand(search_agents) * (upperbound - lowerbound)

    # 计算每个搜索代理的适应度
    fit = np.zeros(search_agents)
    for i in range(search_agents):
        L = X[i, :]
        # print(L)
        # print(type(L))
        fit[i] = fitness_func(L)

    # 存储最佳解
    best_so_far = []

    for t in range(1, max_iterations + 1):
        # 更新最佳候选解
        best_idx = np.argmin(fit)
        Xbest = X[best_idx, :]
        fbest = fit[best_idx]
        best_so_far.append(fbest)

        # 第一阶段：河流或池塘中的河马位置更新（探索）
        for i in range(search_agents // 2):
            Dominant_hippopotamus = Xbest
            I1, I2 = random.randint(1, 3), random.randint(1, 3)
            Ip1 = np.random.randint(0, 2, size=2)
            RandGroup = np.random.choice(search_agents, size=random.randint(1, search_agents), replace=False)

            if len(RandGroup) > 1:
                MeanGroup = np.mean(X[RandGroup], axis=0)
            else:
                MeanGroup = X[RandGroup[0]]

            Alfa = [
                (I2 * np.random.rand(dimension) + (~Ip1[0])).astype(float),
                (2 * np.random.rand(dimension) - 1).astype(float),
                np.random.rand(dimension).astype(float),
                (I1 * np.random.rand(dimension) + (~Ip1[1])).astype(float),
                np.random.rand()
            ]

            A = Alfa[random.randint(0, 4)]
            B = Alfa[random.randint(0, 4)]
            X_P1[i, :] = X[i, :] + random.random() * (Dominant_hippopotamus - I1 * X[i, :])
            T = np.exp(-t / max_iterations)

            if T > 0.6:
                X_P2[i, :] = X[i, :] + A * (Dominant_hippopotamus - I2 * MeanGroup)
            else:
                if random.random() > 0.5:
                    X_P2[i, :] = X[i, :] + B * (MeanGroup - Dominant_hippopotamus)
                else:
                    X_P2[i, :] = np.random.uniform(lowerbound, upperbound)

            X_P2[i, :] = np.clip(X_P2[i, :], lowerbound, upperbound)

            F_P1, F_P2 = fitness_func(X_P1[i, :]), fitness_func(X_P2[i, :])

            if F_P1 < fit[i]:
                X[i, :], fit[i] = X_P1[i, :], F_P1
            if F_P2 < fit[i]:
                X[i, :], fit[i] = X_P2[i, :], F_P2

        # 第二阶段：河马防御捕食者（探索）
        for i in range(search_agents // 2, search_agents):
            predator = np.random.uniform(lowerbound, upperbound,size=(1, dimension))
            predator = predator[0]
            F_HL = fitness_func(predator)
            distance2Leader = np.abs(predator - X[i, :])
            b, c, d, l = np.random.uniform(2, 4), np.random.uniform(1, 1.5), np.random.uniform(2, 3), np.random.uniform(
                -2 * np.pi, 2 * np.pi)
            RL = 0.05 * levy(search_agents, dimension, 1.5)

            if fit[i] > F_HL:
                X_P3[i, :] = RL[i] * predator + (b / (c - d * np.cos(l))) * (1 / distance2Leader)
            else:
                X_P3[i, :] = RL[i] * predator + (b / (c - d * np.cos(l))) * (
                        1 / (2 * distance2Leader + np.random.rand(dimension)))

            X_P3[i, :] = np.clip(X_P3[i, :], lowerbound, upperbound)
            F_P3 = fitness_func(X_P3[i, :])

            if F_P3 < fit[i]:
                X[i, :], fit[i] = X_P3[i, :], F_P3

        # 第三阶段：河马逃离捕食者（开发）
        for i in range(search_agents):
            LO_LOCAL = lowerbound / t
            HI_LOCAL = upperbound / t
            Alfa = [
                (2 * np.random.rand(dimension) - 1).astype(float),
                np.random.rand(),
                np.random.randn(dimension)
            ]
            D = Alfa[random.randint(0, 2)]
            X_P4[i, :] = X[i, :] + random.random() * (LO_LOCAL + D * (HI_LOCAL - LO_LOCAL))
            X_P4[i, :] = np.clip(X_P4[i, :], lowerbound, upperbound)
            F_P4 = fitness_func(X_P4[i, :])

            if F_P4 < fit[i]:
                X[i, :], fit[i] = X_P4[i, :], F_P4

        # 打印迭代信息
        print(f'Iteration {t}: Best Cost = {best_so_far[-1]}')

    Best_score = fbest
    Best_pos = Xbest
    HO_curve = best_so_far

    return Best_score, Best_pos, HO_curve


######################################### end ####################################################

def main():
    Fun_name = 'F4'
    SearchAgents = 16
    Max_iterations = 1000

    # 获取函数信息
    lowerbound, upperbound, dimension, fitness = fun_info(Fun_name)

    # 运行优化算法
    Best_score, Best_pos, HO_curve = HO(SearchAgents, Max_iterations, lowerbound, upperbound, dimension, fitness)

    # 显示结果
    print(f"The best solution obtained by HO for {Fun_name} is : {Best_pos}")
    print(f"The best optimal value of the objective function found by HO for {Fun_name} is : {Best_score}")

    # 绘制结果
    plt.figure()
    plt.semilogy(HO_curve, color='#b28d90', linewidth=2)
    plt.xlabel('Iteration')
    plt.ylabel('Best score obtained so far')
    plt.box(True)
    plt.rcParams['font.family'] = 'Times New Roman'  # 设置全局字体为Times New Roman
    plt.legend(['HO'])
    plt.show()


if __name__ == "__main__":
    main()

如在求函数y=x^2+y^2+10的最小值时，设置的最大迭代次数为1000，dimension即m=2，算法运行结果如下。