冠豪猪优化算法

在优化算法的不断演进中，受自然界生物行为启发的算法逐渐崭露头角，我们之前了解过 蚁群算法、 金豺算法等，接下来介绍一种于2024年M Abdel-Basset等人提出了冠豪猪优化算法（Crested Porcupine Optimizer, CPO），参考文章《 Crested Porcupine Optimizer: A new nature-inspired metaheuristic》，我来讲解以下冠豪猪优化算法的原理及其应用

冠豪猪优化算法原理

冠豪猪从低到高的防御机制依次为视觉、听觉、气味和物理攻击，前两者对应 CPO 的探索行为，后两者对应开发行为。CPO 引入循环种群缩减技术，模拟部分冠豪猪受威胁才激活防御机制的情况，以提升收敛速度和种群多样性。

防御机制

以下四种防御机制分别对应图中的A B C D四个区域

1. 视觉防御：当冠豪猪察觉到危险，它会竖起长而锋利的刺，让自己看起来体型更大，以此吓退捕食者。
2. 声音防御：冠豪猪会发出各种声音，如咆哮、跺脚、磨牙等，随着捕食者靠近，声音强度会增加。
3. 气味防御：若视觉和声音防御未能奏效，冠豪猪会从背部腺体分泌出恶臭的化学气味，以此驱赶捕食者。
4. 物理攻击防御：这是冠豪猪最强大的防御手段，当其他策略都失败时，它会用背上短而粗的刺向后撞击攻击者。

CPO算法的数学模型与实现

初始化阶段

CPO算法在开始时，会在搜索空间内随机初始化一组个体（候选解）。每个个体的维度根据问题的要求在一定范围内随机生成，公式为：

$${\vec{X}}_i=\vec{L}+\vec{r}\times (\vec{U}-\vec{L})|i=1,\mathrm{2...},N$$

其中$N$是个体数量，${\vec{X}}_i$是第$i$个候选解，$\vec{L}$和$\vec{U}$分别是搜索范围的下限和上限，$\vec{r}$是在$0$到$1$之间随机生成的向量。

循环种群缩减技术

CPO算法引入了循环种群缩减技术（Cyclic Population Reduction，CPR），该技术模拟了只有部分冠豪猪在受到威胁时才会激活防御机制的现象。在优化过程中，通过调整种群大小来平衡探索和开发的能力。种群大小的更新公式为：$$N=N_{min}+\left(N^{\prime }-N_{min}\right)\times \left(1-\left(\frac{t\%\frac{T_{max}}{T}}{\frac{T_{max}}{T}}\right)\right)$$ 其中，$T$是循环次数，$t$是当前迭代次数，$T_{max}$是最大迭代次数，$N_{min}$是最小种群大小，$N^{\prime }$是最大种群大小。

取 $T_{max}=120$，$N_{min}=40$，$N^{\prime }=100$，如下图分别为$T=2$、$T=3$、$T=5$的种群数量循环

探索阶段

1. 第一种防御策略：模拟冠豪猪竖起刺的行为，通过生成基于正态分布的随机值来决定捕食者（在算法中可理解为搜索方向）是靠近还是远离冠豪猪（当前解）。当随机值属于$[-1,1]$时，捕食者靠近；否则，远离。$$\overrightarrow{x_i^{t+1}}=\overrightarrow{x_i^t}+{\tau }_1\times \left|2\times {\tau }_2\times \overrightarrow{x_{CP}^t}-\overrightarrow{y_i^t}\right|$$ 其中,$\overrightarrow{x_{CP}^t}$是当前最优解，${\tau }_1$是基于正态分布的随机数，${\tau }_2$是在$[0,1]$区间的随机值，$\overrightarrow{y_i^t}$是代表捕食者位置的向量，其更新公式为$$\overrightarrow{y_i^t}=\frac{\overrightarrow{x_i^t}+\overrightarrow{x_r^t}}{2}$$ 其中，$r$为$[1,N]$的随机数

2. 第二种防御策略：模拟冠豪猪发出声音的行为，通过以下公式更新解的位置：$$\overrightarrow{x_i^{t+1}}=\left(1-\overrightarrow{U_1}\right)\times \overrightarrow{x_i^t}+\overrightarrow{U_1}\times \left(\overrightarrow{y}+{\tau }_3\times \left(\overrightarrow{x_{r1}^r}-\overrightarrow{x_{r2}^t}\right)\right)$$ 其中，$r1$和$r2$是在$[1,N]$之间的随机整数，${\tau }_3$是在$[0,1]$区间的随机值，$\overrightarrow{U_1}$是包含$0$和$1$的随机二进制向量，用于模拟声音对捕食者行为的不同影响。
   此公式模拟了冠豪猪发出的声音对捕食者与自身距离的影响。这种影响分为三种情况：
   \quad (1) 声音较弱时，捕食者继续向冠豪猪移动；
   \quad (2) 声音稍大时，捕食者可能会原地不动；
   \quad (3) 声音很大时，捕食者因恐惧而远离。

开发阶段

1. 第三种防御策略：模拟冠豪猪分泌气味进行防御的行为，通过以下公式更新解的位置：$$\overrightarrow{x_i^{t+1}}=\left(1-\overrightarrow{U_1}\right)\times \overrightarrow{x_i^{\prime }}+\overrightarrow{U_1}\times \left(\overrightarrow{x_1^{\prime }}+S_i^t\times \left(\overrightarrow{x_{r_2}^{\prime }}-\overrightarrow{x_{r_1}^{\prime }}\right)-{\tau }_3\times \overrightarrow{\delta }\times {\gamma }_t\times S_i^t\right)$$其中，涉及多个参数，如$r_3$是$[1,N]$之间的随机整数，$\overrightarrow{\delta }$控制搜索方向，${\gamma }_t$是防御因子，$S_i^t$是气味扩散因子。气味扩散因子的大小决定了开发的范围，小的值使开发更靠近随机选择的解，大的值则使开发远离该解。

   • $r_3$是在$[1,N]$之间的随机数，用于引入随机性，帮助算法在搜索空间中探索不同的区域。

   • $\overrightarrow{\delta }$是用于控制搜索方向的参数$$\overrightarrow{\delta }=\{\begin{array}{ll}+1,& \text{if }\overrightarrow{rand}\le 0.5\\ -1,& \text{Else}\end{array}$$ $\overrightarrow{rand}$是一个在$0$到$1$之间随机生成的数值组成的向量，它决定了$\overrightarrow{\delta }$的值，进而影响搜索方向。

   • $\overrightarrow{x_i^t}$表示第$i$个个体在迭代$t$时的位置，是算法当前正在更新的解的位置。

   • ${\gamma }_t$是防御因子$${\gamma }_t=2\times rand\times (1-\frac{t}{t_{max}}{)}^{\frac{t}{t_{max}}}$$ 其中$rand$是一个在$[0,1]$区间的随机生成的向量，$t$是当前迭代次数，$t_{max}$是最大迭代次数。${\gamma }_t$的值会随着迭代次数的变化而变化，用于调整算法的搜索强度和方向。
     

   • ${\tau }_3$是一个在$[0,1]$区间内的随机值，为算法的更新过程引入更多的随机性。

   • $S_i^t$是 气味扩散因子，由公式 $$S_t^i=\exp \left(\frac{f(x_t^i)}{{\sum }_{k=1}^{N}f(x_t^k)+ϵ}\right)$$ 其中，$f(x_t^i)$代表第$i$个个体在迭代$t$时的目标函数值，$ϵ$是一个很小的值，用于避免分母为零的情况，$N$是种群大小。$S_i^t$的值反映了个体在种群中的相对适应度，其取值范围在$0.3$到$2.6$之间，不同的值对应不同的气味扩散率，从而影响算法的开发行为。

   • $\overrightarrow{U_1}$向量用于模拟该策略中可能出现的三种场景：当$\overrightarrow{U_1}$等于$0$时，意味着冠豪猪停止气味扩散，因为捕食者因恐惧停止移动，此时捕食者与冠豪猪的距离保持不变；当$\overrightarrow{U_1}$等于$1$时，表示冠豪猪会显著地释放气味，因为捕食者已经靠近；当$\overrightarrow{U_1}$是$0$和$1$的组合时，捕食者会与冠豪猪保持安全距离，此时不需要大量释放气味 。

2. 第四种防御策略：模拟冠豪猪的物理攻击行为，当捕食者离得很近时，冠豪猪会用短而厚的刺攻击捕食者。在物理攻击中，两个身体强烈融合，代表一维的非弹性碰撞。通过以下公式更新解的位置：$$\overrightarrow{x_i^{t+1}}=\overrightarrow{x_{CP}^t}+\left(\alpha \left(1-{\tau }_4\right)+{\tau }_4\right)\times \left(\delta \times \overrightarrow{x_{CP}^t}-\overrightarrow{x_i^t}\right)-{\tau }_5\times \delta \times {\gamma }_t\times {\overrightarrow{F}}_i^t$$
   其中，

   • $\overrightarrow{x_{CP}^t}$是当前获得的最佳解，代表冠豪猪的位置;
   • $\overrightarrow{x_i^t}$是第$i$个个体在迭代$t$时的位置，代表捕食者的位置；
   • $\alpha$是收敛速度因子；
   • ${\tau }_4$和${\tau }_5$是在$[0,1]$区间的随机值；
   • $\overrightarrow{F_i^t}$是冠豪猪对第$i$个捕食者的平均作用力 ,由物理中的非弹性碰撞定律计算得出，公式为：
     $$\begin{array}{c}\overrightarrow{F_i^t}=\overrightarrow{{\tau }_6}\times \frac{m_i\times \left(\overrightarrow{v_i^{t+1}}-\overrightarrow{v_i^t}\right)}{\Delta t}\\ m_i=e^{\frac{f\left(\overrightarrow{x_i^t}\right)}{{\sum }_{k=1}^N\left(\overrightarrow{x_k^t}\right)+ϵ}}\\ \overrightarrow{v_i^t}=\overrightarrow{x_i^t}\\ \overrightarrow{v_i^{t+1}}=\overrightarrow{x_r^t}\end{array}$$
     其中，$m_t^i$是第$i$个个体（捕食者）在迭代$t$时的质量，它通过目标函数值$f(\cdot )$来计算，$ϵ$是一个很小的值，用于避免分母为零；$\overrightarrow{v_i^{t+1}}$是第$i$个个体在下一次迭代$t+1$时的最终速度，它是基于从当前种群中随机选择的一个解来赋值的；$\overrightarrow{v_i^t}$是第$i$个个体在迭代$t$时的初始速度，这里直接等于当前个体的位置$\overrightarrow{x_i^t}$；$\Delta t$是当前迭代次数；$\overrightarrow{{\tau }_6}$是在$[0,1]$区间的随机值的组成的向量。
     在原始公式中，平均作用力$\overrightarrow{F_i^t}$的计算是通过将分子除以当前迭代次数$\Delta t$，随着优化过程中迭代次数线性增加，该平均作用力的影响会逐渐减小。如图10所示，这种小的值对CPO算法的性能可能不利，因为它可能无法帮助算法充分探索当前最优解周围的不同区域以找到更好的解。
     
     为解决这一问题，作者对公式进行了更新，删除了分子，仅依赖分母，得到新公式：
     $$\overrightarrow{F_i^t}=\overrightarrow{{\tau }_6}\times m_i\times \left(\overrightarrow{v_i^{t+1}}-\overrightarrow{v_i^t}\right)$$
     这种调整有助于在搜索空间中创建更广泛的值范围，使算法能够更全面地检查当前最优解周围的区域，加快向接近最优解的收敛速度。如图11所示，通过这种方式，算法的开发阶段能够更集中地在当前最优解附近的区域进行搜索，椭圆区域代表了搜索区域在搜索空间内的集中分布，从而提高算法找到更优解的能力。
     

（注：这里有一根很奇怪的点，第三阶段的气味扩散因子$S_i^t$和第四阶段的质量$m_i$，使用的公式是相同的，但是一个使用exp(·)，另一个使用$e^{(\cdot )}$，同时使用的目标函数也都是$f(\cdot )$的形式，没有区分，可能会产生误解）

算法流程

CPO算法在每一代中，根据随机生成的数值决定采用探索阶段还是开发阶段的策略来更新解的位置。更新后的解通过目标函数进行评估，以判断是否达到最优解或满足终止条件。整个算法流程不断迭代，直到达到最大函数评估次数或满足其他终止条件为止。文章中给出伪代码

代码

文章中作者给出了代码来源https://drive.mathworks.com/sharing/24c48ec7-bfd5-4c22-9805-42b7c394c691/
代码中使用了IEEE Congress on Evolutionary Computation (CEC )会议发布的测试函数套件（如 CEC 2014、CEC 2017、CEC 2020、CEC 2022）为优化算法的性能评估提供了标准化的工具。

main.m 主程序

1. 初始化参数，设置种群数量130，最大迭代次数$T_{max}=1000000$

%%
clear all
clc
CP_no=130; % Number of search agents (crested porcupines)
Tmax=1000000; % Maximum number of Function evaluations

2. 设置算法独立运行30次和Fun_id 30组参数选择（需要多次独立运行以评估算法的稳定性和性能）

RUN_NO=30; %% Number of independent runs
Fun_id=[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];

3. 指定测试函数CEC 2017

fhd=str2func('cec17_func'); %%Default Benchmark
cec=2;
if cec==1 %% CEC-2014
    fhd=str2func('cec14_func');
    benchmarksType='cec14_func';
elseif cec==2 %% CEC-2017
    fhd=str2func('cec17_func');
    benchmarksType='cec17_func';
elseif cec==3 %% CEC-2020
    fhd=str2func('cec20_func');
    benchmarksType='cec20_func';
elseif cec==4 %% CEC-2022
    fhd=str2func('cec22_func');
    benchmarksType='cec22_func';
end

4. 针对30组参数，每组运行30次（使用parfor函数多进程加速循环运行），保存Convergence_curve和fitness，并对每组参数绘制30次运行平均收敛曲线

for i=1:30  % 针对30组参数
    tic;
    if cec==2 && i==2
        continue;
    elseif cec==3 && i>10
        return
    elseif cec==4 && i>12
        return
    end
    parfor j=1:RUN_NO  % 每组运行30次
        [lb,ub,dim]=Get_Functions_detailsCEC(Fun_id(i));  % 获取当前参数
        fobj=Fun_id(i);
        [Best_score,Best_pos,Convergence_curve(j,:)]=CPO(CP_no,Tmax,ub,lb,dim,fobj,fhd);  % 使用CPO算法
        fitness(1,j)=Best_score;  % 保存适应度参数
    end
    time=toc;
    fprintf(['benchmark   \t',num2str(cec),'\t','Function_ID\t',num2str(i),'\tAverage Fitness:',num2str(mean(fitness(1,:)),20),'\tAverage Time:',num2str(time/RUN_NO,5),'\n']);
    %% Drawing Convergence Curve %% 绘制收敛曲线
    figure(i)
    h=semilogy(mean(Convergence_curve),'-<','MarkerSize',8,'LineWidth',1.5);
    h.MarkerIndices = 1000:20000:Tmax;
    xlabel('Function Evaluation');
    ylabel('Best Fitness obtained so-far');
    axis tight
    grid off
    box on
    legend({'CPO'});
    clear Convergence_curve;
end

以第一组参数为例，展示收敛曲线结果

Get_Functions_detailsCEC.m 函数

函数包含了论文中表1、表2和表3中基准测试函数的全部信息和实现

• lb 是下界：lb=[lb_1, lb_2, ..., lb_d]
• ub 是上界：ub=[ub_1, ub_2, ..., ub_d]
• dim 是变量的数量（问题的维度）

function [lb,ub,dim] = Get_Functions_detailsCEC(F)

函数中包含30组参数，除第11组外，其余均一样

CPO.m 冠豪猪优化算法主体函数

该函数展示了文章中讲解的CPO算法是如何运行的

1. 函数声明

function [Gb_Fit, Gb_Sol, Conv_curve] = CPO(Pop_size, Tmax, ub, lb, dim, fobj, fhd)

• 输入

  • Pop_size：种群大小
  • Tmax：最大迭代次数
  • ub：变量的上界
  • lb：变量的下界
  • dim：问题的维度
  • fobj：目标函数
  • fhd：测试函数CEC

• 输出

  • Gb_Fit：全局最优解的适应度值
  • Gb_Sol：全局最优解
  • Conv_curve：收敛曲线

2. 定义输出参数的类型和大小

Gb_Fit=zeros(1,dim); % A vector to include the best-so-far solution
Gb_Sol=inf; % A Scalar variable to include the best-so-far score
Conv_curve=zeros(1,Tmax);

3. 给定控制参数

N=Pop_size; %% Is the initial population size. 初始种群数目
N_min=120; %% Is the minimum population size. 最小种群数目
T=2; %% The number of cycles  % 种群数量循环次数
alpha=0.2; %% The convergence rate  % 收敛率
Tf=0.8; %% The percentage of the tradeoff between the third and fourth defense mechanisms 第三种和第四种防御机制之间的权衡比例

4. 初始化冠豪猪的初始位置

X=initialization(Pop_size,dim,ub,lb); % Initialize the positions of crested porcupines
t=0; %% Function evaluation counter 

5. 针对初始位置，求解每一个冠豪猪的对cec2017的适应度，并寻找适应度最小的个体Gb_sol

for i=1:Pop_size
    %% Test suites of CEC-2014, CEC-2017, CEC-2020, and CEC-2022
    fitness(i)=feval(fhd, X(i,:)',fobj);
end
% 寻找最小适应Gb_Fit及其对应的个体Gb_sol
[Gb_Fit,index]=min(fitness);
Gb_Sol=X(index,:);    
%% A new array to store the personal best position for each crested porcupine
Xp=X;
opt=fobj*100; %% Well-known optimal fitness 最优适应度

6. 进行CPO优化求解
   我会按照自己的理解将代码中变量名和文章中的公式符号尽量对应起来

   • 判定条件：迭代次数$t$超出最大迭代次数$T_{max}$ 或者 当前最小适应Gb_Fit等于最优适应度opt，否则就循环进行CPO优化
     

   • 代码中变量 r2 对应 第四步物理防御策略中 ${\tau }_4$，是$[0,1]$之间的随机值，各次迭代间r2值随机互不相同，但一次迭代中针对每一只个体保持r2参数相同
     

   • 针对每一只冠豪猪个体进行CPO优化，更新位置信息，具体步骤如下

     • 通过随机值 rand < rand之间比较，判断当前是 探索阶段 还是 开发阶段

       • 探索阶段——依旧是通过rand < rand判断是 第一防御策略 还是 第二防御策略
         

         • 第一防御策略
           • 首先生成捕食者位置 y 对应符号 $\overrightarrow{y_i^t}$，使用公式$$\overrightarrow{y_i^t}=\frac{\overrightarrow{x_i^t}+\overrightarrow{x_r^t}}{2}$$
           • 而后更新第i个冠豪猪位置 X(i,:) 对应符号 $\overrightarrow{x_i^t}$，使用公式 $$\overrightarrow{x_i^{t+1}}=\overrightarrow{x_i^t}+{\tau }_1\times \left|2\times {\tau }_2\times \overrightarrow{x_{CP}^t}-\overrightarrow{y_i^t}\right|$$ 符号${\tau }_1$和${\tau }_2$分别使用函数randn和rand代替，符号$\overrightarrow{x_{CP}^t}$对应Gb_sol。
             但是代码中的更新方式 顺序依赖性 有种弊端：就是先更新第一只冠豪猪位置，然后使用已经更新的第一冠豪猪位置来更新第二只，而不是使用原来旧的冠豪猪位置来更新
         • 第二防御策略
           • 同样是首先生成捕食者位置 y 对应符号 $\overrightarrow{y_i^t}$
           • 而后更新第i个冠豪猪位置 X(i,:) 对应符号 $\overrightarrow{x_i^t}$，使用公式和文章中不同 $$\overrightarrow{x_i^{t+1}}=\overrightarrow{U_1}\times \overrightarrow{x_i^t}+\left(1-\overrightarrow{U_1}\right)\times \left(\overrightarrow{y}+{\tau }_3\times \left(\overrightarrow{x_{r1}^r}-\overrightarrow{x_{r2}^t}\right)\right)$$ 程序中将$\overrightarrow{U_1}$和$(1-\overrightarrow{U_1})$的位置互换了，$\overrightarrow{U_1}$本身是binary随机向量，所以互换后和原式作用基本相同，丝毫不影响程序运行，但对应文章中$\overrightarrow{U_1}$取值0或1的含义解释可能不同

       • 开发阶段——首先生成捕食者位置 Yt 对应符号 $\overrightarrow{y_i^t}$，而后使用文章中公式(8)生成U2对应符号$\overrightarrow{\delta }$，之后计算S同时代表第三步中${\tau }_3\times \overrightarrow{y_i^t}$和第四步中${\tau }_5\times \overrightarrow{y_i^t}$。
         最后依旧是通过rand < Tf 判断是 第三防御策略 还是 第四防御策略，这里的变量Tf取值为0.8，在给定控制参数步骤中设置
         

         • 第三防御策略
           • 计算St代表符号$S_i^t$气味扩散因子 ，使用公式$$S_t^i=\exp \left(\frac{f(x_t^i)}{{\sum }_{k=1}^{N}f(x_t^k)+ϵ}\right)$$ 此时的$f(\cdot )$是原始设定的适应度函数fintness。
           • 更新更新第i个冠豪猪位置 X(i,:) 对应符号 $\overrightarrow{x_i^t}$，这里代码有错误，直接使用Yt来充当符号${\gamma }_t$，这是一个很显著的错误，在第四防御策略中也出现了这样的错误。出错位置如下
             
         • 第四防御策略
           • 变量Mt代表符号$m_i$
           • 变量vt代表符号$\overrightarrow{v_i^t}$
           • 变量Vtp代表符号$\overrightarrow{v_i^{t+1}}$
           • 变量Ft代表符号$\overrightarrow{F_i^t}$
           • 更新更新第i个冠豪猪位置 X(i,:) 对应符号 $\overrightarrow{x_i^t}$，这里代码有错误，同样也直接使用Yt来充当符号${\gamma }_t$

   • 判断位置更新后是否超出范围[lb,ub]，如若超出则随机生成新位置
     

   • 更新全局和局部最优解，并保存最小适应度fitness作为收敛曲线
     

CPO进行路径规划

绑定资源中有代码，之后有时间再来对路径规划代码详细介绍

参考资料：
[1] Mohamed Abdel-Basset, Reda Mohamed, Mohamed Abouhawwash. Crested Porcupine Optimizer: A new nature-inspired metaheuristic[J]. Knowledge-Based Systems, 2024, 284: 111257.