智能优化算法——哈里斯鹰算法（Matlab代码实现）
目录
1 写在前面
先搞笑一下：
在自然界中，哈里斯鹰会利用其犀利的双眼侦查环境、追踪猎物。但是，在茫茫的亚利桑那州南部地区，有时候日子并不好过。在沙漠地区，其常常需要花费几个小时来等待，观察，并追踪猎物。
哈里斯鹰算法(Harris Hawks Optimizer，HHO)是由Heidari、Mirjalili等人于2019年提出的一种元启发式算法，不得不说Mirjalili教授是真滴厉害！哈里斯之鹰主要生活在美国亚利桑那州南部地区，其之所以与众不同，是因为它会与群体中的其他家庭成员一起进行独特的合作觅食活动，而其他种类的猛禽则通常独自捕食猎物。正因如此，哈里斯鹰独特的群体捕食行为才非常适合被模拟成一种群智能优化过程。
讨论一个算法怎么样，首先得看看他的有点吧：该算法有较强的全局搜索能力，并且需要调节的参数较少的优点。
2 哈里斯鹰优化数学模型
2.1 位置更新公式
哈里斯鹰优化(HHO)是一种近两年被提出的基于种群的优化思想,其核心在于模拟自然界中哈里斯鹰群的动态捕猎模式3]。在HHO的探索阶段，为了模拟哈里斯鹰侦查猎物的状态,鹰群的个体位置按下式随机更新：
其中q、r1、r2、r3、r4均为[0,1]内的随机数，ub、lb为搜索空间上、下限；Xrand 为一随机个体位置；Xrabbit 为猎物位置，Xave为种群内所有个体的平均位置。r3是一个比例系数，一旦r4的值接近1，就会进一步增加该策略的随机性。类似于鲸鱼优化算法，绝对值中的内容可看作两个体间的相对距离；r1为随机的比例系数，其为哈里斯鹰的栖息提供多样化趋势并使其能够探索不同区域的特征空间。
2.2 猎物能量减少的线性递减公式
与其他的群智能算法相似,HHO的探索与开发阶段的转换是通过一个模拟猎物能量减少的线性递减下式控制的。哈里斯鹰可以根据猎物(小兔兔)的逃逸能量来进行不同状态的转换。在猎物的逃脱过程中，其逃逸能量E 将大大降低：
式中,E是猎物的初始能量,每次迭代时在(-1，1)范围内更新; Maxiter表示最大迭代代数。
鉴于不同猎物间的逃逸能量存在差异，故原文令E0(逃逸能量初始值)在算法迭代过程中于[-1,1]内随机变化。关于这个逃逸能量初始值，原文给了这样的解释：
“
When the value of E0 decreases from 0 to -1, the rabbit is physically flagging, whilst when the value of E0 increases from 0 to 1, it means that the rabbit is strengthening. The dynamic escaping energy E has a decreasing trend during the iterations.”
当E0由0减小至-1时(E0<0)，此时兔兔“身体越加疲乏”？其实可以理解为兔兔不断逃跑耗费了很多能量，所以越来越虚；当E0由0增加至1时(E0>0)，兔兔处于恢复能量阶段。当|E |≥1时，哈里斯鹰搜索不同的区域以进一步探索猎物的位置，此时对应全局搜索阶段；当|E |<1时，哈里斯鹰对相邻的解进行局部勘探，故而对应局部开发阶段。
另一方面，HHO中的开发阶段按照哈里斯鹰的不同追逐模式可以分为以下4种策略：
在说四种策略之前，要先提一下攻击的前提。假设r 为兔兔逃脱概率，r <0.5时为成功逃脱，相反r ≥0.5为逃脱失败。通常情况下，哈里斯鹰会以强硬或轻柔的围攻方式来捕获猎物(一柔一刚，不得不佩服作者的想象力)。这围攻方式意味着哈里斯鹰将根据猎物的剩余能量从不同方向轻柔或强硬地攻击猎物。在实际情况下，哈里斯鹰会越来越接近预期的猎物，并通过突袭而增加了合作杀死猎物的机会。一段时间后，猎物将损失越来越多的能量，这时哈里斯鹰便加强围攻过程，从而捕获猎物。在这个过程中，逃逸能量E 的作用不言而喻。原文假设|E |≥0.5时，进行轻柔围攻；当|E |<0.5时，进行强硬围攻。
2.3 软包围
可以通过下式表示：
其中J 为兔兔逃跑过程中的跳跃距离，J=2*(1-rand)。
敌强我退，敌疲我进，老毛的大策略都用上了，赞！
兔兔能量太多，先来软滴：
2.4 硬包围
2.5 采取渐进式快速俯冲的软包围策略
LF(·)为莱维飞行的数学表达式,在已知群的维度D的基础上,更新策略可表达为:
2.6 采取渐进式快速俯冲的硬包围策略
3 matlab代码实现
3.1 代码
%%===欢迎关注公众号：电力系统与算法之美===
clear
close all
clc
SearchAgents_no = 30 ; % 种群规模
dim = 10 ; % 粒子维度
Max_iter = 1000 ; % 迭代次数
ub = 5 ;
lb = -5 ;
%% 初始化猎物位置和逃逸能量
Rabbit_Location=zeros(1,dim);
Rabbit_Energy=inf;
%% 初始化种群的位置
Positions= lb + rand(SearchAgents_no,dim).*(ub-lb) ;
Convergence_curve = zeros(Max_iter,1);
%% 开始循环
for t=1:Max_iter
for i=1:size(Positions,1)
FU=Positions(i,:)>ub;FL=Positions(i,:)<lb;Positions(i,:)=(Positions(i,:).*(~(FU+FL)))+ub.*FU+lb.*FL;
fitness=sum(Positions(i,:).^2);
if fitness<Rabbit_Energy
Rabbit_Energy=fitness;
Rabbit_Location=Positions(i,:);
end
end
E1=2*(1-(t/Max_iter));
%% 鹰群的个体位置位置更新
for i=1:size(Positions,1)
E0=2*rand()-1; %-1<E0<1
Escaping_Energy=E1*(E0);
if abs(Escaping_Energy)>=1
%%
q=rand();
rand_Hawk_index = floor(SearchAgents_no*rand()+1);
X_rand = Positions(rand_Hawk_index, :);
if q<0.5
%%
Positions(i,:)=X_rand-rand()*abs(X_rand-2*rand()*Positions(i,:));
elseif q>=0.5
Positions(i,:)=(Rabbit_Location(1,:)-mean(Positions))-rand()*((ub-lb)*rand+lb);
end
elseif abs(Escaping_Energy)<1
%%
%% phase 1
r=rand();
if r>=0.5 && abs(Escaping_Energy)<0.5
Positions(i,:)=(Rabbit_Location)-Escaping_Energy*abs(Rabbit_Location-Positions(i,:));
end
if r>=0.5 && abs(Escaping_Energy)>=0.5
Jump_strength=2*(1-rand());
Positions(i,:)=(Rabbit_Location-Positions(i,:))-Escaping_Energy*abs(Jump_strength*Rabbit_Location-Positions(i,:));
end
%% phase 2
if r<0.5 && abs(Escaping_Energy)>=0.5
Jump_strength=2*(1-rand());
X1=Rabbit_Location-Escaping_Energy*abs(Jump_strength*Rabbit_Location-Positions(i,:));
if sum(X1.^2)<sum(Positions(i,:).^2)
Positions(i,:)=X1;
else
beta=1.5;
sigma=(gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta);
u=randn(1,dim)*sigma;v=randn(1,dim);step=u./abs(v).^(1/beta);
o1=0.01*step;
X2=Rabbit_Location-Escaping_Energy*abs(Jump_strength*Rabbit_Location-Positions(i,:))+rand(1,dim).*o1;
if (sum(X2.^2)<sum(Positions(i,:).^2))% improved move
Positions(i,:)=X2;
end
end
end
if r<0.5 && abs(Escaping_Energy)<0.5
Jump_strength=2*(1-rand());
X1=Rabbit_Location-Escaping_Energy*abs(Jump_strength*Rabbit_Location-mean(Positions));
if sum(X1.^2)<sum(Positions(i,:).^2)
Positions(i,:)=X1;
else
beta=1.5;
sigma=(gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta);
u=randn(1,dim)*sigma;v=randn(1,dim);step=u./abs(v).^(1/beta);
o2=0.01*step;
X2=Rabbit_Location-Escaping_Energy*abs(Jump_strength*Rabbit_Location-mean(Positions))+rand(1,dim).*o2;
if (sum(X2.^2)<sum(Positions(i,:).^2))% improved move
Positions(i,:)=X2;
end
end
end
%%
end
end
Convergence_curve(t)=Rabbit_Energy;
if mod(t,50)==0
display(['At iteration ', num2str(t), ' the best fitness is ', num2str(Rabbit_Energy)]);
end
end
figure('unit','normalize','Position',[0.3,0.35,0.4,0.35],'color',[1 1 1],'toolbar','none')
subplot(1,2,1);
x = -5:0.1:5;y=x;
L=length(x);
f=zeros(L,L);
for i=1:L
for j=1:L
f(i,j) = x(i)^8+y(j)^8;
end
end
surfc(x,y,f,'LineStyle','none');
xlabel('x_1');
ylabel('x_2');
zlabel('F')
title('Objective space')
subplot(1,2,2);
semilogy(Convergence_curve,'Color','r','linewidth',1.5)
title('Convergence_curve')
xlabel('Iteration');
ylabel('Best score obtained so far');
axis tight
grid on
box on
legend('HHO')
display(['The best solution obtained by HHO is : ', num2str(Rabbit_Location)]);
display(['The best optimal value of the objective funciton found by HHO is : ', num2str(Rabbit_Energy)]);
3.2 结果