RHSSA1.m

function [GBestX,GBestF,curve]=RHSSA1(fhd,Number,Max_iter,lb,ub,dim,func_num)  

% fhd=str2func('cec17_func');
% Number=100;
% dim=30;
% Max_iter=500;
% lb=-100.*ones(1,dim);
% ub=100.*ones(1,dim);
% func_num=8;

curve=zeros(1,Max_iter);
ST = 0.6;
PD = 0.3;
SD = 0.1;

PDNumber = round(Number*PD); 
SDNumber = round(SD*Number);
if(max(size(ub)) == 1)
   ub = ub.*ones(1,dim);
   lb = lb.*ones(1,dim);  
end
%初始化

pop0=initialization(Number,dim,ub,lb);
pop=pop0;
fitness = zeros(1,Number);
for i = 1:Number
   fitness(i) =  feval(fhd,pop0(i,:)',func_num);
end
 [fitness, index]= sort(fitness);%
 for i=1:Number
     pop(i,:)=pop0(index(i),:);
 end
%  b=feval(fhd,pop(1,:)',func_num);
 
GBestF = fitness(1); 
GBestX=pop(1,:);

% sz=50;
LP=40;

Z1=ones(1,Max_iter);%%记录扰动成功与失败的标志——柯西变异
Z2=ones(1,Max_iter);%%记录扰动成功与失败的标志---高斯变异
 
pop_next=zeros(Number,dim);
fit_pop_next=zeros(1,Number);
for t= 1: Max_iter
    
    XBestF = fitness(1);
    XBest=pop(1,:);

    WorstF = fitness(end);
    Xworst=pop(end,:);
    %% 划分发现者和跟随者
    pop_temp1=zeros(Number,dim);
    fit_p_temp1=zeros(1,Number);
    pop_temp1_1=zeros(Number,dim);
    fit_p_temp1_1=zeros(1,Number);
     pop_temp2=zeros(Number,dim);
   
     D=zeros(Number,Number);
     Dmmin=zeros(1,Number);
     div=zeros(1,Number);
     C=zeros(1,Number);
     Gt=zeros(1,dim);
 index1=zeros(1,Number);
     for i=2:Number
        for j=1:i-1
             D(i,j)=sqrt(sum((pop(i,:)-pop(j,:)).^2));
      end
     end
%%记录相对距离     
 for i=2:Number
     A=D(i,1:i-1);
     [a,~]=min(A);
     Dmmin(i)=a;
     
 end
 Dmmin(1)=max(Dmmin);
 
 %% 排序
    [Dmin,qt] =sort(Dmmin,'descend');%%相对距离降序排列
    %%得到适应度值排序为i的个体的相对距离排序为q(i)
    qq=zeros(1,Number);
    for j=1:Number
        i=1;
         while (qt(i)~=j)
            i=i+1;
        end
        qq(j)=i;
    end
        
        
        
    for i=1:Number
    div(i)=qq(i)+i;
    end
    [div,div1]=sort(div);

div1_1=div1(1:PDNumber);
div1_2=div1(PDNumber+1:Number);
%% div1_1前PDxN记录发现者的下标
%%  div1_2记录着成为跟随者的个体下标
%% 因为最后还要是排序，所以用w=1从1开始接收，之后再进行排序
m=1;w=1;
sum1=0;
for m=1:PDNumber
   sum1=fitness(div1_1(m))+sum1;
end
m=1;
%%求出新选择的发现者对应的权重
for m=1:PDNumber
C(div1_1(m))=(fitness(end)-fitness(div1_1(m)))/(PDNumber*fitness(end)-sum1);
end
% sum(C)
%%求出重心
Gt=zeros(1,dim);
for m=1:PDNumber
    Gt=Gt+C(div1_1(m)).*pop(div1_1(m),:);
end
w=1;m=1;
%exp(-(i)/(r1*M))
while(m<=PDNumber)
    %%先进行位置更新，再求出反向解，再保留适应度值较低的个体让让pop_next接收
    R2=rand;
    
    if (R2<ST)
        ro=rand();
    pop_temp1(div1_1(m),:)=pop(div1_1(m),:).*exp(-div1_1(m)/(ro*Max_iter));
    
    else
    pop_temp1(div1_1(m),:)=pop(div1_1(m),:)+randn(1,dim).*ones(1,dim);
%     m=m+1;
    end
    pop_temp1( div1_1(m), : ) = Bounds( pop_temp1(div1_1(m),:), lb, ub );
    fit_p_temp1( div1_1(m)) = feval( fhd,pop_temp1(div1_1(m), : )',func_num);   
    %%求出反向解
    pop_temp1_1(div1_1(m),:)=2*Gt-pop_temp1(div1_1(m),:);
    %%控制边界
    for k=1:dim
        if(pop_temp1_1(div1(m),k)>ub(k))
            pop_temp1_1(div1(m),k)=Gt(1,k)+rand*(ub(k)-Gt(k));
        end
        if(pop_temp1_1(div1(m),k)<lb(k))
            pop_temp1_1(div1(m),k)=lb(k)+rand*(Gt(k)-lb(k));
        end
    end
    fit_p_temp1_1(div1_1(m))=feval( fhd,pop_temp1_1(div1_1(m), : )',func_num);
    
    %%保留适应度较好的个体
    if(feval(fhd,pop_temp1(div1_1(m),:)',func_num)<feval(fhd,pop_temp1_1(div1_1(m),:)',func_num))
        pop_next(w,:)=pop_temp1(div1_1(m),:);
    else
        pop_next(w,:)=pop_temp1_1(div1_1(m),:);
    end
     fit_pop_next(w,:)=feval( fhd,pop_next(div1_1(m), : )',func_num);
    w=w+1;
    m=m+1;
    %% 求出加权重心    
end
    %% 进行跟随者位置更新
    for m=PDNumber+1:Number
        %%因为采用m,所以此处用总的div1
        if (m>Number/2)
            pop_next(w,:)=randn(1,dim).*exp(pop(end,:)-pop(div1(m),:));
             pop_next(w, : ) = Bounds( pop_next(w,:), lb, ub );
        else
            A = ones(1,dim);
          for a = 1:dim
            if(rand()>0.5)
                A(a) = -1;
            end
          end 
          AA = A'*inv(A*A');  
          %%因为适应度值最小的发现者相对距离设最大，所以此处可直接用pop(1,:)是一样的结果
          pop_next(w,:)= pop(1,:) + abs(pop(div1(m),:) - pop(1,:)).*AA';
          pop_next(w, : ) = Bounds( pop_next(w,:), lb, ub );
        end
        fit_pop_next(w,:)=feval(fhd,pop_next(w, : )',func_num);
        w=w+1;%w及时加1，以便于接收解
    end
    
    
    for i=1:Number
      if(feval(fhd,pop_next(i,:)',func_num)<GBestF)
          GBestX=pop_next(i,:);
          GBestF=feval(fhd,pop_next(i,:)',func_num);
      end
    end
  
%% 警戒值位置更新
  Temp = randperm(Number);
   SDchooseIndex = Temp(1:SDNumber); 
   for i = 1:SDNumber
       if(fitness(SDchooseIndex(i))>XBestF)
           pop_next(SDchooseIndex(i),:) = pop(1,:) + randn().*abs(pop(SDchooseIndex(i),:) - pop(1,:));
          
       elseif(fitness(SDchooseIndex(i))== XBestF)
           K = 2*rand() -1;
           pop_next(SDchooseIndex(i),:) = pop(SDchooseIndex(i),:) + K.*(abs(pop(SDchooseIndex(i),:) - pop(end,:))./(fitness(SDchooseIndex(i)) - fitness(end) + 10^-8));
          
       end
         pop_next(SDchooseIndex(i), : ) = Bounds( pop_next(SDchooseIndex(i),:) ,lb, ub);
         fit_pop_next(SDchooseIndex(i),:)=feval( fhd,pop_next(SDchooseIndex(i), : )',func_num);
   end
    for i=1:Number
      if(feval(fhd,pop_next(i,:)',func_num)<GBestF)
          GBestX=pop_next(i,:);
          GBestF=feval(fhd,pop_next(i,:)',func_num);
      end
    end
   
   XBest_Temp1=XBest+randn(1,dim).*XBest;%高斯变异
   XBest_Temp1 = Bounds( XBest_Temp1, lb, ub );
   r = tan((rand(1,dim) - 0.5)*pi);%柯西随机数
    XBest_Temp2=XBest+r.*XBest;
   XBest_Temp2 = Bounds( XBest_Temp2, lb, ub );
    if(feval(fhd,XBest_Temp1',func_num)<feval(fhd,XBest_Temp2',func_num))
        XBest_Final1=XBest_Temp1;
        fit_XBest_Final1=feval(fhd,XBest_Temp1',func_num);
        Z2(1,t)=0;
    else
        XBest_Final1=XBest_Temp2;
        fit_XBest_Final1=feval(fhd,XBest_Temp2',func_num);
        Z1(1,t)=0;
    end
    if t<=LP
        p=1/2;
    else
    p=sum(Z1(t-LP:t-1))/(sum(Z1(t-LP:t-1))+sum(Z2(t-LP:t-1)));
    end
    R1=rand();
    if p<=R1
        XBest_Final= XBest_Temp1;
         fit_XBest_Final=feval(fhd,XBest_Temp1',func_num);
    else
         XBest_Final= XBest_Temp2;
         fit_XBest_Final=feval(fhd,XBest_Temp2',func_num);
    end
    if(fit_XBest_Final<GBestF)
        GBestF=fit_XBest_Final;
        GBestX=XBest_Final;
    end
    if(fit_XBest_Final<XBestF)
        XBestF=fit_XBest_Final;
        XBest=XBest_Final;
    end
             

        

  %%更新全局最优解
  for i=1:Number
      if(feval(fhd,pop_next(i,:)',func_num)<GBestF)
          GBestX=pop_next(i,:);
          GBestF=feval(fhd,pop_next(i,:)',func_num);
      end
  end
  %赋值排序，准备下一次迭代。
   pop0=pop_next;
   for i=1:Number
      fitness(i)=feval(fhd,pop0(i,:)',func_num);
   end
  [fitness,index2]=sort(fitness);
   for i=1:Number
      pop(i,:)=pop0(index2(i),:);
   end
%   a=feval(fhd,pop(1,:)',func_num)
   curve(t)=GBestF;
%    pause(10);
%    t
%    GBestF
 end
end