冠豪猪优化算法(CPO)（附开源完整MATLAB代码）2024最新算法！全网最全面讲解！

1.算法简介

引用自：Crested Porcupine Optimizer: A new nature-inspired metaheuristic.

冠豪猪优化器(CPO)是由Abdel-Basset等人在2024年初提出，并首次发表在中科院1区的SCI顶级期刊《Knowledge-Based Systems》上，是目前最新也是最具竞争力的优化算法之一。其基本思想是模仿冠豪猪在寻找食物、避免天敌以及在复杂环境中生存的策略。不同于现有的智能优化算法，CPO模拟防御机制的灵感在所有现有元启发式技术中被认为是独一无二的。

2.核心思想与优势

流程简介

• 防御策略模拟：CPO算法模仿豪猪的四种防御策略，包括视觉恐吓、声音恐吓、气味攻击和物理攻击，分别对应于算法的探索和剥削阶段，这种基于防御机制的算法在元启发式领域被认为是独一无二的。
• 动态交互与更新：豪猪之间进行动态交互，包括信息共享、竞争和协作，以促进种群向更优解进化。
• 群体动态调整策略：CPO引入了群体缩减策略CPR，通过在优化过程中让一些CP离开群体，随后再次加入群体以提高多样性。加速收敛速度的同时也提升了群体的多样性。
• 改进率较高：CPO可以实现比竞争对手优化器更高的改进率，CEC2017的百分比高达83%，CEC2017的百分比为70%，CEC2020的百分比为90%，六个实际工程问题的百分比为100%。

3.算法主要步骤

3.1 初始化阶段

优化问题由多个维度(d)构成，每个维度均需通过优化技术进行最小化或最大化处理。CPO算法通过在每个问题的搜索空间内随机初始化它们的维度来优化这些问题，具体使用以下公式进行初始化：

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

其中，$N^{\prime }$表示个体数量，${\vec{X}}_i$是搜索空间中的第i个候选解，$\vec{L}$和$\vec{U}$分别是搜索范围的下界和上界，$\vec{r}$是在0和1之间随机初始化的向量。

3.2 循环群体缩减技术

CPR是一种新颖的方法，旨在保持群体多样性的同时加速收敛速度。该策略通过在优化过程中让一些CP离开群体，以加速收敛速度，随后再次加入群体以提高多样性，避免陷入局部最小值。循环群体大小减少的数学模型如下：

$$N=N_{min}+(N^{\prime }-N_{min})\times \left(1-{\left(\frac{t\%T_{max}}{T}\right)}^{\frac{T}{T_{max}}}\right)$$

其中，T是确定循环次数的变量，t是当前函数评估，$T_{max}$是最大函数评估次数，%表示取余运算符，$N_{min}$是新生成群体中个体的最小数量。

循环种群缩减技术下的种群规模描述如下图所示，:

在图(a)中，当T设定为2时，优化过程开始时所有CP都被激活。随后的过程如下：

1. 群体大小随着当前函数评估的增加而逐渐减少，直至达到40，这代表了第一个循环的探索与缩减阶段。
2. 然后，群体大小再次恢复到最大，标志着第一个循环的结束。
3. 接下来，第二个循环开始，群体大小再次经历从最大到最小的渐变过程，直至优化过程结束。
4. 因为T设置为2，这意味着整个优化过程将经历两次这样的循环。

通过这种方式，群体大小呈现出先最大化，然后逐渐减少至$N_{\text{min}}$的周期性变化。这一周期性变化在优化过程中重复出现，直到完成了T次循环，如图所示。

3.3 探索阶段

3.3.1 第一个防御策略

当CP意识到捕食者时，它开始抬起并扇动其刺毛，给人一种更大的印象。这种行为的数学模拟如下：

$${\vec{x}}_i^{t+1}={\vec{x}}_i^t+{\tau }_1\times \sqrt{2}\times {\tau }_2\times ({\vec{x}}_{CP}^t-{\vec{y}}_i^t)$$

其中，${\vec{x}}_i^t$是第t次迭代中第i个CP的位置，${\vec{x}}_{CP}^t$是当前最佳CP的位置，${\tau }_1$和${\tau }_2$是随机数，${\vec{y}}_i^t$是通过以下公式计算的：

$${\vec{y}}_i^t=\frac{{\vec{x}}_i^t+{\vec{x}}_k^t}{2}$$

这里，${\vec{x}}_k^t$是随机选择的另一个CP的位置。

模拟第一防御机制的搜索历史如下图所示：

该图模拟了该防御机制的搜索历史，可以看出，它为CPO提供了尽可能覆盖整个搜索空间的强大潜力，从而避免陷入局部极小值。

3.3.2 第二个防御策略

CP使用声音方法发出噪音并威胁捕食者。随着捕食者的接近，豪猪的声音变得越来越响亮。为了数学上模拟这种行为，提出了以下公式：

$${\vec{x}}_i^{t+1}=(1-{\vec{U}}_1)\times {\vec{x}}_i^t+{\vec{U}}_1\times (\vec{y}+{\tau }_3\times ({\vec{x}}_{r1}^t-{\vec{x}}_{r2}^t))$$

其中，${\vec{U}}_1$是一个0到1之间的随机向量，${\tau }_3$是一个随机数，${\vec{x}}_{r1}^t$和${\vec{x}}_{r2}^t$是随机选择的两个CP的位置。

第二种防御机制中捕食者和CP的行为模拟如下图所示：

3.4 开发阶段

3.4.1 第三个防御策略

在这种策略中，CP分泌恶臭气味，扩散在其周围区域，防止捕食者接近。为了数学上模拟这种行为，提出了以下公式：

$${\vec{x}}_i^{t+1}=(1-{\vec{U}}_1)\times {\vec{x}}_i^t+{\vec{U}}_1\times ({\vec{x}}_{r1}^t+S_t^i\times ({\vec{x}}_{r2}^t-{\vec{x}}_{r3}^t)-{\tau }_3\times \vec{\delta }\times {\gamma }_t\times S_t^i)$$

其中，$S_t^i$是一个适应度函数，$\vec{\delta }$是一个随机向量，${\gamma }_t$是一个时间相关因子，${\vec{x}}_{r1}^t$、${\vec{x}}_{r2}^t$和${\vec{x}}_{r3}^t$是随机选择的三个CP的位置。

气味扩散系数的各种值与第三种防御机制的搜索历史仿真图如下图所示：
在这个策略中的另一个重要术语是气味扩散因子（St），它负责在优化过程中控制气味的扩散速率。根据上式，这个因子的取值范围在0.3到2.6之间，如左图所示。如果这个因子的值很小，那么扩散速率将非常低，因此搜索的开发过程将集中在随机选择的解附近，这在右图中用椭圆形状高亮表示。相反，如果这些值较高，那么气味的扩散速率会很高，从而搜索的开发过程会在远离随机选定解的位置进行，如右图中非椭圆高亮区域所示。

3.4.2 第四个防御策略

最后一个策略是身体攻击。在数学公式中表达其身体攻击行为如下：

$${\vec{x}}_i^{t+1}={\vec{x}}_{CP}^t+(\alpha (1-{\tau }_4)+{\tau }_4)\times (\delta \times {\vec{x}}_{CP}^t-{\vec{x}}_i^t)-{\tau }_5\times \delta \times {\gamma }_t\times {\vec{F}}_t^i$$

其中，$\alpha$是一个控制参数，${\tau }_4$和${\tau }_5$是随机数，$\delta$是一个随机向量，${\vec{F}}_t^i$是一个基于适应度的因子。

CPO的时间复杂度分析表明，其性能主要取决于问题的维度和最大函数评估次数，这与其他元启发式算法一致。总体而言，CPO在大O记号下的时间复杂度为$O(T_{max}d)$。

仿真第四种防御机制的搜索历史如下图所示，其中椭圆形表示开发区域在搜索空间内的小区域上的色散。

4.伪代码与核心代码

伪代码与核心代码

4.1 CPO的伪代码

1. Set parameters $N$, $T_{max}$, $T_f$, $T$, $N_{min}$;
2. Initialize the solutions’ positions randomly, $x_i$, $i=1,2,...,N$;
3. While $t<T_{max}$
   1. Evaluate fitness values for the candidate solutions;
   2. Determine the best solution so far, $x_{CP}^t$;
   3. Update $\gamma$ using Eq. (9).
   4. Update the population size using Eq. (3)
   5. For $i=1:N$:
      1. Update $m$, $S.F.$, $\overline{\delta }$;
      2. Generate two random numbers, ${\tau }_8$ and ${\tau }_9$
      3. If ${\tau }_8<{\tau }_9$: % Exploration phase
         1. Generate two random numbers, ${\tau }_6$ and ${\tau }_7$
         2. If ${\tau }_6<{\tau }_7$: % First defense mechanism
            • Apply Eq.(4)
         3. Else: % Second defense mechanism
            • Apply Eq.(6)
      4. Else: % Exploitation phase
         1. Generate a random number, ${\tau }_{10}$
         2. If ${\tau }_{10}<T_f$: % Third defense mechanism
            • Apply Eq.(7)
         3. Else: % Fourth defense mechanism
            • Apply Eq.(11)
      5. End If
      6. If $f(x_i^{t+1})>f(x_i^t)$
         • $x_i^{t+1}=x_i^t$
      7. End If
   6. t = t + 1
   7. End For
4. End While
5. Return $x_{CP}^t$

4.2 CPO的核心代码

原论文中给出了一份论文作者写的开源代码，链接在这: 作者开源代码
但是这份代码主函数使用并行加速，容易报错，为此在博客也开源了一份小编修改后的完整代码，可以流畅运行，完全免费，可放心自取。
链接放在这里：
完整代码链接：https://pan.baidu.com/s/10sDMGkaeerzLWp9qH47kaA
提取码：9999

function [Score, Best_pos, CPO_curve] = CPO(Search_Agents, Max_iterations, Lowerbound, Upperbound, dimensions, objective)
% CPO: Crested Porcupine Optimization 算法
% 输入:
%   Search_Agents: 搜索代理数量
%   Max_iterations: 最大迭代次数
%   Lowerbound: 搜索空间下界
%   Upperbound: 搜索空间上界
%   dimensions: 问题维度
%   objective: 目标函数句柄

% 初始化搜索边界
Lowerbound = ones(1, dimensions) .* Lowerbound;
Upperbound = ones(1, dimensions) .* Upperbound;

% 初始化控制参数
N_min = 120; % 种群规模最小值
T = 2; % 循环数量
alpha = 0.2; % 收敛率
Tf = 0.8; % 第三和第四防御机制之间的权衡百分比

% 初始化种群和适应度
X = initialization(Search_Agents, dimensions, Upperbound, Lowerbound);
fitness = zeros(1, Search_Agents);
CPO_curve = zeros(1, Max_iterations);

% 计算初始适应度
for i = 1:Search_Agents
    fitness(i) = objective(X(i,:));
end

% 初始化全局最优
[Score, index] = min(fitness);
Best_pos = X(index, :);

% 存储每个峰冠豪猪个人最佳位置
Xp = X;

% 主循环
t = 0;
while t < Max_iterations
    r2 = rand;
    for i = 1:Search_Agents
        U1 = rand(1, dimensions) > rand;
        
        if rand < rand % 探索阶段
            if rand < rand % 第一防御机制
                rand_index = randi(Search_Agents);
                y = (X(i,:) + X(rand_index,:)) / 2;
                X(i,:) = X(i,:) + (randn) .* abs(2*rand*Best_pos - y);
            else % 第二防御机制
                rand_index1 = randi(Search_Agents);
                rand_index2 = randi(Search_Agents);
                y = (X(i,:) + X(rand_index1,:)) / 2;
                X(i,:) = (U1) .* X(i,:) + (1-U1) .* (y + rand*(X(rand_index1,:) - X(rand_index2,:)));
            end
        else
            Yt = 2*rand*(1-t/Max_iterations)^(t/Max_iterations);
            U2 = rand(1,dimensions) < 0.5*2-1;
            S = rand*U2;
            if rand < Tf % 第三防御机制
                St = exp(fitness(i)/(sum(fitness)+eps));
                S = S.*Yt.*St;
                rand_index1 = randi(Search_Agents);
                rand_index2 = randi(Search_Agents);
                rand_index3 = randi(Search_Agents);
                X(i,:) = (1-U1).*X(i,:) + U1.*(X(rand_index1,:) + St*(X(rand_index2,:) - X(rand_index3,:)) - S);
            else % 第四防御机制
                Mt = exp(fitness(i)/(sum(fitness)+eps));
                vt = X(i,:);
                rand_index = randi(Search_Agents);
                Vtp = X(rand_index,:);
                Ft = rand(1,dimensions).*(Mt*(-vt+Vtp));
                S = S.*Yt.*Ft;
                X(i,:) = (Best_pos + (alpha*(1-r2)+r2)*(U2.*Best_pos-X(i,:))) - S;
            end
        end
        
        % 边界检查
        X(i,:) = max(X(i,:), Lowerbound);
        X(i,:) = min(X(i,:), Upperbound);
        
        % 计算新适应度
        new_fitness = objective(X(i,:));
        
        % 更新个体最优和全局最优
        if new_fitness < fitness(i)
            Xp(i,:) = X(i,:);
            fitness(i) = new_fitness;
            if new_fitness < Score
                Best_pos = X(i,:);
                Score = new_fitness;
            end
        else
            X(i,:) = Xp(i,:);
        end
    end
    
    t = t + 1;
    CPO_curve(t) = Score;
    
    % 动态调整种群规模
    New_Search_Agents = fix(N_min + (Search_Agents - N_min) * (1 - (rem(t, Max_iterations/T) / (Max_iterations/T))));
    
    if New_Search_Agents < Search_Agents
        [~, sorted_indexes] = sort(fitness);
        X = X(sorted_indexes(1:New_Search_Agents), :);
        Xp = Xp(sorted_indexes(1:New_Search_Agents), :);
        fitness = fitness(sorted_indexes(1:New_Search_Agents));
        Search_Agents = New_Search_Agents;
    end
end

end

% 辅助函数：初始化种群
function X = initialization(SearchAgents_no, dim, ub, lb)
    Boundary_no = size(ub, 2);
    if Boundary_no == 1
        X = rand(SearchAgents_no, dim) .* (ub - lb) + lb;
    else
        for i = 1:dim
            ub_i = ub(i);
            lb_i = lb(i);
            X(:, i) = rand(SearchAgents_no, 1) .* (ub_i - lb_i) + lb_i;
        end
    end
end

5.基准函数测试

论文中对算法在多个基准函数上进行了测试，测试结果为：

完整代码（免费自取）

链接：https://pan.baidu.com/s/10sDMGkaeerzLWp9qH47kaA
提取码：9999

本篇文章不靠卖代码赚取收益，麻烦给个点赞和关注，后续还会有开源的免费优化算法及其代码，栓Q！同时如果大家有想要的算法可以在评论区打出，如果有空的话我可以帮忙复现