野马优化算法WHO Matlab代码

算法主体

% Max_iter: 最大迭代, N: populatoin size, Convergence_curve:收敛曲线
function [Convergence_curve,gBest,gBestScore]=WHO(N,Max_iter,lb,ub,dim,fobj)

 if size(ub,2)==1
     ub=ones(1,dim)*ub;
     lb=ones(1,dim)*lb;
 end
 PS=0.2;     % Stallions Percentage 
 PC=0.13;    % Crossover Percentage
 NStallion=ceil(PS*N); % number Stallion头领个数
 Nfoal=N-NStallion; % 马驹个数

Convergence_curve = zeros(1,Max_iter);
% gBest=zeros(1,dim);
% gBestScore=inf;  % inf为无穷大
 
% create initial population
empty.pos=[]; % empty是一个结构体，empty包含pos，pos：[]
empty.cost=[];

% group包含所有马驹
group=repmat(empty,Nfoal,1); % 生成 Nfoal 行，1列，元素为empty的堆叠矩阵

% 给group的foal马驹赋值
for i=1:Nfoal 
   group(i).pos=lb+rand(1,dim).*(ub-lb);
   group(i).cost=fobj(group(i).pos);
end

% 给头领赋值
Stallion=repmat(empty,NStallion,1);

for i=1:NStallion 
   Stallion(i).pos=lb+rand(1,dim).*(ub-lb);
   Stallion(i).cost=fobj(Stallion(i).pos);
 
end

ngroup=length(group); % group的个数
a=randperm(ngroup); % 生成ngroup个元素的整数行向量，元素从1到ngroup随机排列
group=group(a); % 打乱group的原排序

% 分组
i=0;
k=1;
for j=1:ngroup
  i=i+1;    
  Stallion(i).group(k)=group(j); % 给每一个group随机赋值
  if i==NStallion
    i=0;
    k=k+1;
  end
end

Stallion=exchange(Stallion); % 变换首领
[~,index]=min([Stallion.cost]);

WH=Stallion(index); % 全局最小
   gBest=WH.pos;
   gBestScore=WH.cost;

Convergence_curve(1)=WH.cost; % 收敛曲线
l=2; % 循环计数器，全局变量
while l<Max_iter+1
    TDR=1-l*((1)/Max_iter);

    for i=1:NStallion
    
        ngroup=length(Stallion(i).group); % 组的成员数
        [~,index]=sort([Stallion(i).group.cost]); % index为从小到大的序号
        Stallion(i).group=Stallion(i).group(index);
        % 以上两行可以这样缩写
        %[Stallion(i).group,index]=sort([Stallion(i).group.cost]);
   
        for j=1:ngroup
    
            if rand>PC
                z=rand(1,dim)<TDR; % 逻辑值数组 0，1
                r1=rand;
                r2=rand(1,dim);
                idx=(z==0);
                r3=r1.*idx+r2.*~idx; % ~idx为非idx
                rr=-2+4*r3;
                % 成员都在头领身边
                Stallion(i).group(j).pos= 2*r3.*cos(2*pi*rr).*(Stallion(i).pos-Stallion(i).group(j).pos)+(Stallion(i).pos);
            else
                A=randperm(NStallion);
                A(A==i)=[]; % A=除去A(i)的其它所有元素
                a=A(1);
                c=A(2);

                x1=Stallion(c).group(end).pos;
                x2=Stallion(a).group(end).pos;

                y1=(x1+x2)/2;   % Crossover

                Stallion(i).group(j).pos=y1; % 第i组的第j个成员替换为y1
            end
   

            Stallion(i).group(j).pos=min(Stallion(i).group(j).pos,ub);
            Stallion(i).group(j).pos=max(Stallion(i).group(j).pos,lb);
    
            Stallion(i).group(j).cost=fobj(Stallion(i).group(j).pos);
    

    
        end
    
% end
% 
% for i=1:NStallion

        R=rand;
%     z=rand(1,dim)<TDR;
%             r1=rand;
%             r2=rand(1,dim);
%             idx=(z==0);
%             r3=r1.*idx+r2.*~idx;
%             rr=-2+4*r3;

        if R<0.5
            k= 2*r3.*cos(2*pi*rr).*(WH.pos-(Stallion(i).pos))+WH.pos;
        else
            k= 2*r3.*cos(2*pi*rr).*(WH.pos-(Stallion(i).pos))-WH.pos;
        end

        k=min(k,ub);
        k=max(k,lb);
        fk=fobj(k);
        if fk<Stallion(i).cost
            Stallion(i).pos  =k;
            Stallion(i).cost=fk;
        end
    end
    Stallion=exchange(Stallion);
    [value,index]=min([Stallion.cost]);
   if value<WH.cost
       WH=Stallion(index);
   end
   gBest=WH.pos;
   gBestScore=WH.cost;
   Convergence_curve(l)=WH.cost;
   l = l + 1;
   
end

初始化函数

% 初始化矩阵，上下限为up与down，随机数
function [X]=initialization(N,dim,up,down)

if size(up,1)==1 % 如果up的行数为1
    X=rand(N,dim).*(up-down)+down;
end
if size(up,1)>1
    for i=1:dim
        high=up(i);low=down(i);
        X(:,i)=rand(1,N).*(high-low)+low;
    end
end
end

变换首领函数

%% 更换头领
function Stallion=exchange(Stallion)

nStallion=length(Stallion);



for i=1:nStallion
    
    % 找到一个组的最小值
   [value,index]=min([Stallion(i).group.cost]);
   

   % 更换头领
   if value<Stallion(i).cost
       
       bestgroup=Stallion(i).group(index);
       
       Stallion(i).group(index).pos=Stallion(i).pos;
       Stallion(i).group(index).cost=Stallion(i).cost;
       
       
       Stallion(i).pos=bestgroup.pos;
       Stallion(i).cost=bestgroup.cost;
       
       
       
       
       
   end
    
    
    
    
end

测试函数

function [lb,ub,dim,fobj] = Get_Functions_details(F)


switch F
    case 'F1'              
        fobj = @F1;
        lb=-50;
        ub=50;
        dim=10;
       
        
    case 'F2'            
        fobj = @F2;
        lb=-10;
        ub=10;
        dim=30;

    case 'F3'
        fobj = @F3;
        lb=-65.536;
        ub=65.536;
        dim=2;
        
    case 'F4'
        fobj = @F4;
        lb=-100;
        ub=100;
        dim=30;
        
    case 'F5'
        fobj = @F5;
        lb=-30;
        ub=30;
        dim=30;
        
    case 'F6'
        fobj = @F6;
        lb=-5;
        ub=20;
        dim=2;     
    case 'F7'       %加了约束项的，罚函数
        fobj = @F7;
        lb=-100;
        ub=100;
        dim=30;
     case 'F8'       %加了约束项的，罚函数
        fobj = @F8;

        lb=[0,0,0];
        ub=[1,1,1];
        dim=3;
    
        case 'F9'       %加了约束项的，罚函数
        fobj = @F9;

        lb=[0,0,0];
        ub=2.*[1,1,1];
        dim=3;
end
end

% F1

function o = F1(x)
o=sum(-x.*sin(sqrt(abs(x))));
end

% F2

function o = F2(x)
o=sum(abs(x))+prod(abs(x));
end


% F3

function o = F3(x)
aS=[-32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32 -32 -16 0 16 32;,...
-32 -32 -32 -32 -32 -16 -16 -16 -16 -16 0 0 0 0 0 16 16 16 16 16 32 32 32 32 32];

for j=1:25
    bS(j)=sum((x'-aS(:,j)).^6);
end
o=(1/500+sum(1./([1:25]+bS))).^(-1);
end

% F4

function o = F4(x)
o=max(abs(x));
end

% F5

function o = F5(x)
dim=size(x,2);
o=sum(100*(x(2:dim)-(x(1:dim-1).^2)).^2+(x(1:dim-1)-1).^2);

end

% F6

function o = F6(x)
o=(x(2)-(x(1)^2)*5.1/(4*(pi^2))+5/pi*x(1)-6)^2+10*(1-1/(8*pi))*cos(x(1))+10;
end

function o = F7(x)
o=sum(abs((x+.5)).^2);
%加个约束，罚函数，
alpha=10000;
%不等式约束
g1=(x(1)+x(2)^2-x(3)>0)*(x(1)+x(2)^2-x(3))^2*alpha;    
g2=(x(6)^3+x(7)^2+x(8)^3>0)*(x(6)^3+x(7)^2+x(8)^3)^2*alpha;
o=o+g1+g2;
end

function o = F8(x)
    o1=-(x(1)^2-x(2)^2+x(2)*x(3));    %原始目标函数值

    %加个约束，罚函数，不满足约束进行惩罚
    alpha=10000;
    %不等式约束  要化成小于等于
    c(1)=2*x(1)+x(2)+3*x(3)-6;
    c(2)=x(1)^2+x(1)*x(2)+x(3)*x(2)-x(2)-6;
    %等式约束 
    ceq=0;
%  ceq(1)=x(1)+x(3)-1;

    o=o1+alpha*sum((c>0).*c.^2)+alpha*sum(ceq.^2);
end

function o = F9(x)
    o1=-x(1)-2*x(2)+3*x(3);    %原始目标函数值

    %加个约束，罚函数，不满足约束进行惩罚
    alpha=10000;
    %不等式约束  要化成小于等于
    c(1)=-x(1)-x(2)+3;
    c(2)=-x(3)-x(2)+3;
    %等式约束 
    ceq(1)=x(1)+x(3)-4;
%     ceq(1)=x(1)+x(3)-1;

    o=o1+alpha*sum((c>0).*c.^2)+alpha*sum(ceq.^2);
end

main函数

% Developed in MATLAB R2017b
% Source codes demo version 1.0
% _____________________________________________________

%  Author, inventor and programmer: Iraj Naruei and Farshid Keynia,
%  e-Mail: irajnaruei@iauk.ac.ir , irajnaruei@yahoo.com
% _____________________________________________________
%  Co-author and Advisor: Farshid Keynia
%
%         e-Mail: fkeynia@gmail.com
% _____________________________________________________
%  Co-authors: Amir Sabbagh Molahoseyni
%        
%         e-Mail: sabbagh@iauk.ac.ir
% _____________________________________________________

% You can find the WHO code at 
% _____________________________________________________
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear  
close all
clc

N=30; % Number of search agents

Function_name='F6'; % % Name of the test function that can be from F1 to F13 and cec01 to cec10

Max_iter=500; % Maximum number of iterations

% Load details of the selected benchmark function
[dim,fobj,ub, lb]  = Select_Functions(Function_name);
[Convergence_curve,gBest,gBestScore]=WHO(N,Max_iter,lb,ub,dim,fobj);

display(['The best location of WHO is: ', num2str(gBest)]);
display(['The best fitness of WHO is: ', num2str(gBestScore)]);
plot(Convergence_curve)
gBestScore