%% 智能优化算法——哈里鹰算法（Matlab实现）
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)
        % Check boundries
        FU=Positions(i,:)>ub;FL=Positions(i,:)<lb;Positions(i,:)=(Positions(i,:).*(~(FU+FL)))+ub.*FU+lb.*FL;
        % fitness of locations
        fitness=sum(Positions(i,:).^2);
        % Update the location of Rabbit
        if fitness<Rabbit_Energy
            Rabbit_Energy=fitness;
            Rabbit_Location=Positions(i,:);
        end
    end
    E1=2*(1-(t/Max_iter)); % factor to show the decreaing energy of rabbit
    %% Update the location of Harris' hawks
    for i=1:size(Positions,1)
        E0=2*rand()-1; %-1<E0<1
        Escaping_Energy=E1*(E0);  % escaping energy of rabbit
        if abs(Escaping_Energy)>=1
            %% Exploration:
            % Harris' hawks perch randomly based on 2 strategy:
            q=rand();
            rand_Hawk_index = floor(SearchAgents_no*rand()+1);
            X_rand = Positions(rand_Hawk_index, :);
            if q<0.5
                % perch based on other family members
                Positions(i,:)=X_rand-rand()*abs(X_rand-2*rand()*Positions(i,:));
            elseif q>=0.5
                % perch on a random tall tree (random site inside group's home range)
                Positions(i,:)=(Rabbit_Location(1,:)-mean(Positions))-rand()*((ub-lb)*rand+lb);
            end
        elseif abs(Escaping_Energy)<1
            %% Exploitation:
            % Attacking the rabbit using 4 strategies regarding the behavior of the rabbit
            %% phase 1: surprise pounce (seven kills)
            % surprise pounce (seven kills): multiple, short rapid dives by different hawks
            r=rand(); % probablity of each event
            if r>=0.5 && abs(Escaping_Energy)<0.5 % Hard besiege
                Positions(i,:)=(Rabbit_Location)-Escaping_Energy*abs(Rabbit_Location-Positions(i,:));
            end
            if r>=0.5 && abs(Escaping_Energy)>=0.5  % Soft besiege
                Jump_strength=2*(1-rand()); % random jump strength of the rabbit
                Positions(i,:)=(Rabbit_Location-Positions(i,:))-Escaping_Energy*abs(Jump_strength*Rabbit_Location-Positions(i,:));
            end
            %% phase 2: performing team rapid dives (leapfrog movements)
            if r<0.5 && abs(Escaping_Energy)>=0.5% Soft besiege % rabbit try to escape by many zigzag deceptive motions
                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) % improved move
                    Positions(i,:)=X1;
                else % hawks perform levy-based short rapid dives around the rabbit
                    %Levy flight
                    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% Hard besiege % rabbit try to escape by many zigzag deceptive motions
                % hawks try to decrease their average location with the rabbit
                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) % improved move
                    Positions(i,:)=X1;
                else % Perform levy-based short rapid dives around the rabbit
                    %Levy flight
                    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;
    % Print the progress every 100 iterations
    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)^2+y(j)^2;
    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)]);