【智能优化算法-蝠鲼优化算法】基于蝠鲼优化算法求解多目标优化问题附matlab代码

 1 内容介绍

蝠鲼觅食优化器 (MRFO) 已显示出处理单目标现实世界问题的良好能力,这使其在解决多目标问题中的应用成为一个有趣的方向。因此,本文研究了 MRFO 优化器,以开发一种新的算法来处理多目标工程设计问题。为了实现这一目标,采用精英概念,通过将外部存档集成到标准 MRFO 中来保存 Pareto 解决方案集。该档案也被认为是一个存储库,根据其密度程度从中选择搜索代理以控制蝠鲼种群的收敛性和多样性。我们的算法的效率首先通过对十个测试函数的广泛实验进行验证,在几乎所有情况下,结果在收敛性和多样性方面都非常令人满意。然后,将其应用于四个多目标工程问题,并在解决具有多目标的现实世界问题方面显示出良好的前景。​

2 仿真代码

function G=CreateHypercubes(costs,ngrid,alpha)

    nobj=size(costs,1);
    
    empty_grid.Lower=[];
    empty_grid.Upper=[];
    G=repmat(empty_grid,nobj,1);
    
    for j=1:nobj
        
        min_cj=min(costs(j,:));
        max_cj=max(costs(j,:));
        
        dcj=alpha*(max_cj-min_cj);
        
        min_cj=min_cj-dcj;
        max_cj=max_cj+dcj;
        
        gx=linspace(min_cj,max_cj,ngrid-1);
        
        G(j).Lower=[-inf gx];
        G(j).Upper=[gx inf];
        
    end

end

%___________________________________________________________________%
%  MOMRFO: Multi-Objective Manta Ray Foraging Optimizer for         %
%  Handling engineering design problems                             %
%                                                                   %
%   Main paper:                                                     %
%                                                                   %
%    A. Got, D. Zouache, A. Moussaoui, MOMRFO: Multi-Objective      % 
%    Manta Ray Foraging Optimizer for handling engineering design   %
%    problems, Knowledge Based-Systems, in press,                   %
%    in press, DOI: https://doi.org/10.1016/j.knosys.2021.107880    %       %

format short;
global pop_size ub lb M D R
%% Select the test functions
fun='ZDT1';  
mo_Xrange;
nVar=D;
%% Parameter settings
pop_size=100;       % Population size
MaxIt=2000;         % Maximum Number of Iterations
Ar_size=100;        % The maximum size of archive
alpha=0.1;          % Grid Inflation Parameter
nGrid=30;           % Number of Grids per each Dimension
beta=4;             % Leader Selection Pressure Parameter
gamma=2;            % Extra (to be deleted) archive Member Selection Pressure
%%  The number of independent runs
Nbrun=1; 
%%  START  THE  EXECUTION  OF  THE  ALGORITHM 
for ru=1:Nbrun
%% Initialize the population of Manta Rays  
MantaRay=CreateEmptyMantaRay(pop_size);
for i=1:pop_size
    for j=1:nVar
        MantaRay(i).Position(1,j)=unifrnd(lb(j),ub(j),1);
    end
    MantaRay(i).Cost=mo_test_function(MantaRay(i).Position,fun);  %% fitness function 
    MantaRay(i).Best.Position=MantaRay(i).Position;
    MantaRay(i).Best.Cost=MantaRay(i).Cost;
end
%% Prepare the archive for the first iteration
MantaRay=DetermineDomination(MantaRay);
Archive=GetNonDominatedMantaRay(MantaRay);
Archive_costs=GetCosts(Archive);
G=CreateHypercubes(Archive_costs,nGrid,alpha);
for i=1:numel(Archive)
    [Archive(i).GridIndex Archive(i).GridSubIndex]=GetGridIndex(Archive(i),G);
end
costs=GetCosts(MantaRay);
Archive_costs=GetCosts(Archive);
%% The main loop of MOMRFO algorithm
for t=1:MaxIt  
Coef=t/MaxIt; 
   Leader=SelectLeader(Archive,beta);   % Select Gbesr
   if rand<0.5
      r1=rand;                         
      Beta2=2*exp(r1*((MaxIt-t+1)/MaxIt))*(sin(2*pi*r1));    
      if Coef>rand                                                      
         MantaRay(1).Position=Leader.Position+rand(1,D).*(Leader.Position-MantaRay(1).Position)+Beta2*(Leader.Position-MantaRay(1).Position); %Equation (4)
      else
         IndivRand=rand(1,D).*(ub-lb)+lb;                                
         MantaRay(1).Position=IndivRand+rand(1,D).*(IndivRand-MantaRay(1).Position)+Beta2*(IndivRand-MantaRay(1).Position); %Equation (7)         
      end              
    else 
         Alpha2=2*rand(1,D).*(-log(rand(1,D))).^0.5;           
         MantaRay(1).Position=MantaRay(1).Position+rand(1,D).*(Leader.Position-MantaRay(1).Position)+Alpha2.*(Leader.Position-MantaRay(1).Position); %Equation (1)
   end
   for i=2:pop_size
      Leader=SelectLeader(Archive,beta);   %% Select Gbest
      if rand<0.5
         r1=rand;                         
         Beta2=2*exp(r1*((MaxIt-t+1)/MaxIt))*(sin(2*pi*r1));    
           if Coef>rand                                                      
              MantaRay(i).Position=Leader.Position+rand(1,D).*(MantaRay(i-1).Position-MantaRay(i-1).Position)+Beta2*(Leader.Position-MantaRay(i).Position); %Equation (4)
           else
              IndivRand=rand(1,D).*(ub-lb)+lb;                                
              MantaRay(i).Position=IndivRand+rand(1,D).*(MantaRay(i-1).Position-MantaRay(i).Position)+Beta2*(IndivRand-MantaRay(i).Position);  %Equation (7)       
           end              
      else
          Alpha2=2*rand(1,D).*(-log(rand(1,D))).^0.5;
          MantaRay(i).Position=MantaRay(i).Position+rand(1,D).*(MantaRay(i-1).Position-MantaRay(i).Position)+Alpha2.*(Leader.Position-MantaRay(i).Position); %Equation (1)
      end         
   end
   %% Adjust boundaries if necessary
   for i=1:pop_size
      MantaRay(i).Position=min(max(MantaRay(i).Position,lb),ub);
   end
      S=2;
      for i=1:pop_size   
         Leader=SelectLeader(Archive,beta);   %% Select Gbest
         MantaRay(i).Position=MantaRay(i).Position+S*(rand*Leader.Position-rand*MantaRay(i).Position); %Equation (8)
         MantaRay(i).Position=min(max(MantaRay(i).Position,lb),ub);
         MantaRay(i).Cost=mo_test_function(MantaRay(i).Position,fun);
      end
%% ubdate the archive 
MantaRay=DetermineDomination(MantaRay);
Archive=[Archive;MantaRay];
Archive=DetermineDomination(Archive);
Archive=GetNonDominatedMantaRay(Archive);
for i=1:numel(Archive)
        [Archive(i).GridIndex Archive(i).GridSubIndex]=GetGridIndex(Archive(i),G);
end 
%% Remove EXTRA-solutions from the high crowded hybercubes in the archive
if numel(Archive)>Ar_size
   EXTRA=numel(Archive)-Ar_size;
   Archive=DeleteFromRep(Archive,EXTRA,gamma);   
end
  Archive_costs=GetCosts(Archive);
  G=CreateHypercubes(Archive_costs,nGrid,alpha); 
  costs=GetCosts(MantaRay);
  Archive_costs=GetCosts(Archive);
%% Front represent the final Pareto front 
Front=Archive_costs';
disp(['Iteration ' num2str(t) '  Probleme test ' fun ' Run ' num2str(ru) ': Number of Pareto Solutions = ' num2str(numel(Front(:,1)))]);
pause(0:0.1);
figure(1);
plot(Front(:,1),Front(:,2),'*')
end
end
filename=strcat('MOMRFO_',fun);
save(filename);
disp('END OF EXECUTION, THANKS!');

function obj=mo_test_function(x,fun)
global M
global D

% MOP test functions
% SCH FON  POL KUR
% DT1 ZDT2 ZDT3 ZDT4 ZDT6 
% DTLZ1  DTLZ2  DTLZ3 DTLZ4 DTLZ5 DTLZ6 
% F1 F2 F3 F4  F5 F6 F7 F8 F9
 
 D=length(x);

% SCH
if strcmp(fun,'SCH')
    obj(1) = x(1)^2;
    obj(2) = (x(1)-2)^2;
end
 
% FON
if  strcmp(fun,'FON')
    D=3;
    y1=(x-1/sqrt(3)).^2;
    g1=sum(y1);
    y2=(x+1/sqrt(3)).^2;
    g2=sum(y2);
    obj(1) = 1-exp(-g1);
    obj(2) = 1-exp(-g2);  
end

% POL
if  strcmp(fun,'POL')
   A1=0.5*sin(1)-2*cos(1)+sin(2)-1.5*cos(2);
   A2=1.5*sin(1)-cos(1)+2*sin(2)-0.5*cos(2);
   B1=0.5*sin(x(1))-2*cos(x(1))+sin(x(2))-1.5*cos(x(2));
   B2=1.5*sin(x(1))-cos(x(1))+2*sin(x(2))-0.5*cos(x(2));
   obj(1) = 1+(A1-B1)^2+(A2-B2)^2;
   obj(2) = (x(1)+3)^2+(x(2)+1)^2; 
end

% KUR
if  strcmp(fun,'KUR') 
    x1=x(1:end-1);
    x2=x(2:end);
    y1=-10*exp(-0.2*sqrt(x1.^2+x2.^2));
    obj(1) =sum(y1);
    y2=abs(x).^0.8+5*sin(x.^3);
    obj(2) =sum(y2);   
end

% ZDT1
if strcmp(fun,'ZDT1')
    y=x(2:end);
    gx=1+9*sum(y)/(D-1);
    obj(1)=x(1);
    obj(2)=gx*(1-sqrt(x(1)/gx));
end 

% ZDT2
if strcmp(fun,'ZDT2')
    y=x(2:end);
    gx=1+9*sum(y)/(D-1);
    obj(1)=x(1);
    obj(2)=gx*(1-(x(1)/gx)^2);
  
end

% ZDT3
if strcmp(fun,'ZDT3')
    y=x(2:end);
    gx=1+9*sum(y)/(D-1);
    obj(1)=x(1);
    obj(2)=gx*(1-sqrt(x(1)/gx)-x(1)/gx*sin(10*pi*x(1)));
  
end

% ZDT4
if strcmp(fun,'ZDT4') 
    y=x(2:end);
    gx=1+10*(D-1)+sum(y.*y-10*cos(4*pi*y));
    obj(1) = x(1);
    obj(2) = gx*(1-sqrt(x(1)/gx));
end


% ZDT6
if strcmp(fun,'ZDT6')
    y=x(2:end);
    gx=1+9*(sum(y)/(D-1))^0.25;
    obj(1)=1-exp(-4*x(1))*(sin(6*pi*x(1)))^6;
    obj(2)=gx*(1-(obj(1)/gx)^2);
 
end
 
 %DTLZ1
if  strcmp(fun,'DTLZ1')    
     switch M
         case 2
         xg=x(2:end);
         gx=100*(D-M+1+sum((xg-0.5).^2-cos(20*pi*(xg-0.5))));
         obj(1)=0.5*x(1)*(1+gx);
         obj(2)=0.5*(1-x(1))*(1+gx);     
         case 3
         xg=x(3:end);
         gx=100*(D-M+1+sum((xg-0.5).^2-cos(20*pi*(xg-0.5))));
         obj(1)=0.5*x(1)*x(2)*(1+gx);
         obj(2)=0.5*x(1)*(1-x(2))*(1+gx);
         obj(3)=0.5*(1-x(1))*(1+gx);
     end
end

%DTLZ2
if  strcmp(fun,'DTLZ2') 
     switch M
         case 2
             xg=x(2:end);
             gx=sum((xg-0.5).^2);
             obj(1)=(1+gx)*cos(x(1)*0.5*pi);
             obj(2)=(1+gx)*sin(x(1)*0.5*pi);  
         case 3
             xg=x(3:end);
             gx=sum((xg-0.5).^2);
             obj(1)=(1+gx)*cos(x(1)*0.5*pi)*cos(x(2)*0.5*pi);
             obj(2)=(1+gx)*cos(x(1)*0.5*pi)*sin(x(2)*0.5*pi);
             obj(3)=(1+gx)*sin(x(1)*0.5*pi);
     end
end

 %DTLZ3
if  strcmp(fun,'DTLZ3') 
    switch M
         case 2
             xg=x(2:end);
             gx=100*(D-M+1+sum((xg-0.5).^2-cos(20*pi*(xg-0.5))));
             obj(1)=(1+gx)*cos(x(1)*0.5*pi);
             obj(2)=(1+gx)*sin(x(1)*0.5*pi);
         case 3
             xg=x(3:end);
             gx=100*(D-M+1+sum((xg-0.5).^2-cos(20*pi*(xg-0.5))));
             obj(1)=(1+gx)*cos(x(1)*0.5*pi)*cos(x(2)*0.5*pi);
             obj(2)=(1+gx)*cos(x(1)*0.5*pi)*sin(x(2)*0.5*pi);
             obj(3)=(1+gx)*sin(x(1)*0.5*pi);
    end
end


%DTLZ4
if  strcmp(fun,'DTLZ4') 
      switch M
         case 2
             xg=x(2:end);
             gx=sum((xg-0.5).^2);
             obj(1)=(1+gx)*cos((x(1)^100)*0.5*pi);
             obj(2)=(1+gx)*sin((x(1)^100)*0.5*pi);
         case 3
             xg=x(3:end);
             gx=sum((xg-0.5).^2);
             obj(1)=(1+gx)*cos((x(1)^100)*0.5*pi)*cos((x(2)^100)*0.5*pi);
             obj(2)=(1+gx)*cos((x(1)^100)*0.5*pi)*sin((x(2)^100)*0.5*pi);
             obj(3)=(1+gx)*sin((x(1)^100)*0.5*pi);
      end
end
 
%DTLZ 5
if  strcmp(fun,'DTLZ5') 
     xg=x(3:end);
     gx=sum((xg-0.5).^2);
     Q1 = x(1);
%      Q2=(pi/(4*(1+gx)))*(1+2*gx*x(2)); 
     Q2= 0.5*(1+2*gx*x(2))/(1+gx);
     obj(1)=(1+gx)*cos(Q1*0.5*pi).*cos(Q2*0.5*pi); 
     obj(2)=(1+gx)*cos(Q1*0.5*pi).*sin(Q2*0.5*pi);
     obj(3)=(1+gx)*sin(Q1*0.5*pi); 
end

% DTLZ6
if  strcmp(fun,'DTLZ6') 
     xg=x(3:end);
     gx=sum(xg.^0.1);
     Q1 = x(1);
     %      Q2=(pi/(4*(1+gx)))*(1+2*gx*x(2)); 
     Q2= 0.5*(1+2*gx*x(2))/(1+gx);    
     obj(1)=(1+gx)*cos(Q1*0.5*pi).*cos(Q2*0.5*pi); 
     obj(2)=(1+gx)*cos(Q1*0.5*pi).*sin(Q2*0.5*pi);
     obj(3)=(1+gx)*sin(Q1*0.5*pi); 
end

% DTLZ7
if  strcmp(fun,'DTLZ7') 
     xg=x(3:end);
     gx=1+(9/(D-M+1))*sum(xg);  
     obj(1)=x(1);
     obj(2)=x(2);
     hf=3-(obj(1)/(1+gx))*(1+sin(3*pi*obj(1)))-(obj(2)/(1+gx))*(1+sin(3*pi*obj(2)));
     obj(3)=(1+gx)*hf; 
end

% F1
if  strcmp(fun,'F1') 
    for j = 1:D
        y(j) = x(j)-x(1)^(0.5*(1+3*(j-2)/(D-2)));
    end
    obj(1)=x(1);  obj(2)=1-sqrt(x(1));
    for k= 1:14
        obj(1) = obj(1)+2/14*y(2*k+1)^2;
    end
    for k= 1:15
        obj(2) = obj(2)+2/15*y(2*k)^2;
    end
end
   
% F2
if  strcmp(fun,'F2') 
    for j = 1:D
        y(j) = x(j)-sin(6*pi*x(1)+j*pi/D);
    end
    obj(1)=x(1);  obj(2)=1-sqrt(x(1));
    for k= 1:14
        obj(1) = obj(1)+2/14*y(2*k+1)^2;
    end
    for k= 1:15
        obj(2) = obj(2)+2/15*y(2*k)^2;
    end
end

% F3
if  strcmp(fun,'F3') 
    for j = 1:D
        if mod(j,2)==1
            y(j) = x(j)-0.8*x(1)*cos(6*pi*x(1)+j*pi/D);     
        else
            y(j) = x(j)-0.8*x(1)*sin(6*pi*x(1)+j*pi/D);
        end
    end
    obj(1)=x(1);  obj(2)=1-sqrt(x(1));
    for k= 1:14
        obj(1) = obj(1)+2/14*y(2*k+1)^2;
    end
    for k= 1:15
        obj(2) = obj(2)+2/15*y(2*k)^2;
    end
end

% F4
if  strcmp(fun,'F4') 
    for j = 1:D
        if mod(j,2)==1
            y(j) = x(j)-0.8*x(1)*cos((6*pi*x(1)+j*pi/D)/3);     
        else
            y(j) = x(j)-0.8*x(1)*sin(6*pi*x(1)+j*pi/D);
        end
    end
    obj(1)=x(1);  obj(2)=1-sqrt(x(1));
    for k= 1:14
        obj(1) = obj(1)+2/14*y(2*k+1)^2;
    end
    for k= 1:15
        obj(2) = obj(2)+2/15*y(2*k)^2;
    end
end

% F5
if  strcmp(fun,'F5') 
    for j = 1:D
        if mod(j,2)==1
            y(j) = x(j)-(0.3*x(1)^2*cos(24*pi*x(1)+4*j*pi/D)+0.6*x(1))*cos(6*pi*x(1)+j*pi/D); 
        else
            y(j) = x(j)-(0.3*x(1)^2*cos(24*pi*x(1)+4*j*pi/D)+0.6*x(1))*sin(6*pi*x(1)+j*pi/D); 
        end
    end
    obj(1)=x(1);  obj(2)=1-sqrt(x(1));
    for k= 1:14
        obj(1) = obj(1)+2/14*y(2*k+1)^2;
    end
    for k= 1:15
        obj(2) = obj(2)+2/15*y(2*k)^2;
    end
end

% F6
if  strcmp(fun,'F6') 
    for j = 1:D
         y(j) = x(j)-2*x(2)*sin(2*pi*x(1)+j*pi/D);
    end 
    obj(1)=cos(0.5*x(1)*pi)*cos(0.5*x(2)*pi);  
    obj(2)=cos(0.5*x(1)*pi)*sin(0.5*x(2)*pi); 
    obj(3)=sin(0.5*x(1)*pi);
    for k= 1:3
        obj(1) = obj(1)+2/3*y(3*k+1)^2;
    end
    for k= 1:2
        obj(2) = obj(2)+2/2*y(3*k+2)^2;
    end
    for k= 1:3
        obj(3) = obj(3)+2/3*y(3*k)^2;
    end
end

% F7
if  strcmp(fun,'F7') 
    for j = 1:D
        y(j) = x(j)-x(1)^(0.5*(1+3*(j-2)/(D-2)));
    end
    for j = 1:D
        z(j) = 4*y(j)^2-cos(8*y(j)*pi)+1.0 ;
    end 
    obj(1)=x(1);  obj(2)=1-sqrt(x(1));
    for k= 1:4
        obj(1) = obj(1)+2/4*z(2*k+1);
    end
    for k= 1:5
        obj(2) = obj(2)+2/5*z(2*k);
    end
end

% F8
if  strcmp(fun,'F8') 
    for j = 1:D
        y(j) = x(j)-x(1)^(0.5*(1+3*(j-2)/(D-2)));
    end
    for j = 1:D
        z(j) = cos((20*y(j)*pi)/sqrt(j)) ;
    end 
    a1 = y(3)^2+y(5)^2+y(7)^2+y(9)^2;
    a2 = y(2)^2+y(4)^2+y(6)^2+y(8)^2+y(10)^2;
    b1 = z(3)*z(5)*z(7)*z(9);
    b2 = z(2)*z(4)*z(6)*z(8)*z(10);
    obj(1) = x(1)+2/4*(4*a1-2*b1+2);
    obj(2) = 1-sqrt(x(1))+2/5*(4*a2-2*b2+2); 
end

% F9
if  strcmp(fun,'F9') 
    for j = 1:D
        y(j) = x(j)-sin(6*pi*x(1)+j*pi/D);
    end
    obj(1)=x(1);  obj(2)=1-x(1)^2;
    for k= 1:14
        obj(1) = obj(1)+2/14*y(2*k+1)^2;
    end
    for k= 1:15
        obj(2) = obj(2)+2/15*y(2*k)^2;
    end
end
if  strcmp(fun,'MaF1') 
    
 
            g  = sum((x(3:end)-0.5).^2,2);
            obj = repmat(1+g,1,M) - repmat(1+g,1,M).*fliplr(cumprod([ones(size(g,1),1),x(1:M-1)],2)).*[ones(size(g,1),1),1-x(M-1:-1:1)];
    
end

if  strcmp(fun,'MaF3') 
    
 
           g= 100*(D-M+1+sum((x(M:end)-0.5).^2-cos(20.*pi.*(x(M:end)-0.5)),2));
            obj = repmat(1+g,1,M).*fliplr(cumprod([ones(size(g,1),1),cos(x(1:M-1)*pi/2)],2)).*[ones(size(g,1),1),sin(x(M-1:-1:1)*pi/2)];
            obj = [x(1:M-1).^4,obj(:,M).^2];
end

if  strcmp(fun,'UF1') 
            J1 = 3 : 2 : D;
            J2 = 2 : 2 : D;
            Y  = x - sin(6*pi*repmat(x(:,1),1,D)+repmat(1:D,size(x,1),1)*pi/D);
            obj(:,1) = x(:,1)         + 2*mean(Y(:,J1).^2,2);
            obj(:,2) = 1-sqrt(x(:,1)) + 2*mean(Y(:,J2).^2,2);   
end

% UF2
if  strcmp(fun,'UF2') 
 J1 = 3 : 2 : D;
            J2 = 2 : 2 : D;
            Y       = zeros(size(x));
            X1      = repmat(x(:,1),1,length(J1));
            Y(:,J1) = x(:,J1)-(0.3*X1.^2.*cos(24*pi*X1+4*repmat(J1,size(x,1),1)*pi/D)+0.6*X1).*cos(6*pi*X1+repmat(J1,size(x,1),1)*pi/D);
            X1      = repmat(x(:,1),1,length(J2));
            Y(:,J2) = x(:,J2)-(0.3*X1.^2.*cos(24*pi*X1+4*repmat(J2,size(x,1),1)*pi/D)+0.6*X1).*sin(6*pi*X1+repmat(J2,size(x,1),1)*pi/D); 
            obj(:,1) = x(:,1)         + 2*mean(Y(:,J1).^2,2);
            obj(:,2) = 1-sqrt(x(:,1)) + 2*mean(Y(:,J2).^2,2);  
end

% UF3
if  strcmp(fun,'UF3') 
            J1 = 3 : 2 : D;
            J2 = 2 : 2 : D;
            Y  = x - repmat(x(:,1),1,D).^(0.5*(1+3*(repmat(1:D,size(x,1),1)-2)/(D-2)));                  
            obj(1) = x(:,1)         + 2/length(J1)*(4*sum(Y(:,J1).^2,2)-2*prod(cos(20*Y(:,J1)*pi./sqrt(repmat(J1,size(x,1),1))),2)+2);
            obj(2) = 1-sqrt(x(:,1)) + 2/length(J2)*(4*sum(Y(:,J2).^2,2)-2*prod(cos(20*Y(:,J2)*pi./sqrt(repmat(J2,size(x,1),1))),2)+2);
        
end

% UF4
if  strcmp(fun,'UF4') 
   J1 = 3 : 2 : D;
            J2 = 2 : 2 : D;
            Y  =x - sin(6*pi*repmat(x(:,1),1,D)+repmat(1:D,size(x,1),1)*pi/D);
            hY = abs(Y)./(1+exp(2*abs(Y)));
            obj(1) = x(:,1)      + 2*mean(hY(:,J1),2);
            obj(2) = 1-x(:,1).^2 + 2*mean(hY(:,J2),2);
end

% UF5
if  strcmp(fun,'UF5') 
 J1 = 3 : 2 : D;
            J2 = 2 : 2 : D;
            Y  = x - sin(6*pi*repmat(x(:,1),1,D)+repmat(1:D,size(x,1),1)*pi/D);
            hY = 2*Y.^2 - cos(4*pi*Y) + 1;
            obj(1) = x(:,1)   + (1/20+0.1)*abs(sin(20*pi*x(:,1)))+2*mean(hY(:,J1),2);
            obj(2) = 1-x(:,1) + (1/20+0.1)*abs(sin(20*pi*x(:,1)))+2*mean(hY(:,J2),2);
         
end

% UF6
if  strcmp(fun,'UF6') 
   J1 = 3 : 2 : D;
            J2 = 2 : 2 : D;
            Y  =x - sin(6*pi*repmat(x(:,1),1,D)+repmat(1:D,size(x,1),1)*pi/D);
            obj(1) = x(:,1)   + max(0,2*(1/4+0.1)*sin(4*pi*x(:,1)))+2/length(J1)*(4*sum(Y(:,J1).^2,2)-2*prod(cos(20*Y(:,J1)*pi./sqrt(repmat(J1,size(x,1),1))),2)+2);
            obj(2) = 1-x(:,1) + max(0,2*(1/4+0.1)*sin(4*pi*x(:,1)))+2/length(J2)*(4*sum(Y(:,J2).^2,2)-2*prod(cos(20*Y(:,J2)*pi./sqrt(repmat(J2,size(x,1),1))),2)+2);
       
end

% UF7
if  strcmp(fun,'UF7') 
  J1 = 3 : 2 : D;
            J2 = 2 : 2 : D;
            Y  = x - sin(6*pi*repmat(x(:,1),1,D)+repmat(1:D,size(x,1),1)*pi/D);
            obj(1) = x(:,1).^0.2   + 2*mean(Y(:,J1).^2,2);
            obj(2) = 1-x(:,1).^0.2 + 2*mean(Y(:,J2).^2,2); 
end

% UF8
if  strcmp(fun,'UF8') 
J1 = 4 : 3 : D;
            J2 = 5 : 3 : D;
            J3 = 3 : 3 : D;
            Y  =x - 2*repmat(x(:,2),1,D).*sin(2*pi*repmat(x(:,1),1,D)+repmat(1:D,size(x,1),1)*pi/D);
            obj(1) = cos(0.5*x(:,1)*pi).*cos(0.5*x(:,2)*pi) + 2*mean(Y(:,J1).^2,2);
            obj(2) = cos(0.5*x(:,1)*pi).*sin(0.5*x(:,2)*pi) + 2*mean(Y(:,J2).^2,2);
            obj(3) = sin(0.5*x(:,1)*pi)+ 2*mean(Y(:,J3).^2,2);
          
end

% UF9
if  strcmp(fun,'UF9') 
  J1 = 4 : 3 : D;
            J2 = 5 : 3 : D;
            J3 = 3 : 3 : D;
            Y  = x - 2*repmat(x(:,2),1,D).*sin(2*pi*repmat(x(:,1),1,D)+repmat(1:D,size(x,1),1)*pi/D);
            obj(1) = 0.5*(max(0,1.1*(1-4*(2*x(:,1)-1).^2))+2*x(:,1)).*x(:,2)   + 2*mean(Y(:,J1).^2,2);
            obj(2) = 0.5*(max(0,1.1*(1-4*(2*x(:,1)-1).^2))-2*x(:,1)+2).*x(:,2) + 2*mean(Y(:,J2).^2,2);
            obj(3) = 1-x(:,2)+ 2*mean(Y(:,J3).^2,2);
         
end

% UF10
if  strcmp(fun,'UF10') 
 J1 = 4 : 3 : D;
            J2 = 5 : 3 : D;
            J3 = 3 : 3 : D;
            Y  = x - 2*repmat(x(:,2),1,D).*sin(2*pi*repmat(x(:,1),1,D)+repmat(1:D,size(x,1),1)*pi/D);
            Y  = 4*Y.^2 - cos(8*pi*Y) + 1;
            obj(1) = cos(0.5*x(:,1)*pi).*cos(0.5*x(:,2)*pi) + 2*mean(Y(:,J1),2);
            obj(2) = cos(0.5*x(:,1)*pi).*sin(0.5*x(:,2)*pi) + 2*mean(Y(:,J2),2);
            obj(3) = sin(0.5*x(:,1)*pi)+ 2*mean(Y(:,J3),2);
          
end

% Welded beam design problem
if strcmp(fun,'Weld')
obj(1)=1.10471*x(1)^2*x(2)+0.04811*x(3)*x(4)*(14.0+x(2));
obj(2)=65856000/(30*10^6*x(4)*x(3)^3);
Z=0;
% Penalty constant
lam=10^15;
Q=6000*(14+x(2)/2);
d=sqrt(x(2)^2/4+(x(1)+x(3))^2/4);
J=2*(x(1)*x(2)*sqrt(2)*(x(2)^2/12+(x(1)+x(3))^2/4));
alpha=6000/(sqrt(2)*x(1)*x(2));
beta=Q*d/J;
tau=sqrt(alpha^2+2*alpha*beta*x(2)/(2*d)+beta^2);
sigma=504000/(x(4)*x(3)^2);
tmpf=4.013*(30*10^6)/196;
P=tmpf*sqrt(x(3)^2*x(4)^6/36)*(1-x(3)*sqrt(30/48)/28);
g(1)=tau-13600;
g(2)=sigma-30000;
g(3)=x(1)-x(4);
g(4)=6000-P;
% No equality constraint in this problem, so empty;
geq=[];

% Apply inequality constraints
for k=1:length(g),
    Z=Z+ lam*g(k)^2*getH(g(k));
end
% Apply equality constraints
for k=1:length(geq),
   Z=Z+lam*geq(k)^2*getHeq(geq(k));
end
obj=obj+Z;
end

% Speed reducer design problem
if strcmp(fun,'Speed')
obj(1)=0.7854*x(1)*x(2)^2*(3.3333*x(3)^2+14.9334*x(3)-43.0934)-1.508*x(1)*(x(6)^2+x(7)^2)+7.4777*(x(6)^3+x(7)^3)+0.7854*(x(4)*x(6)^2+x(5)*x(7)^2);
obj(2) =sqrt(((745*x(4))/(x(2)*x(3)))^2+1.69e7)/(0.1*x(6)^3);
Z=0;
% Penalty constant
lam=10^15;
g(1)=27/(x(1)*x(2)^2*x(3))-1;
g(2)=397.5/(x(1)*x(2)^2*x(3)^2)-1;
g(3)=(1.93*x(4)^3)/(x(2)*x(3)*x(6)^4)-1;
g(4)=(1.93*x(5)^3)/(x(2)*x(3)*x(7)^4)-1;
g(5)=((sqrt(((745*x(4))/(x(2)*x(3)))^2+16.9e6))/(110*x(6)^3))-1;
g(6)=((sqrt(((745*x(5))/(x(2)*x(3)))^2+157.5e6))/(85*x(7)^3))-1;
g(7)=((x(2)*x(3))/40)-1;
g(8)=(5*x(2)/x(1))-1;
g(9)=(x(1)/12*x(2))-1;
g(10)=((1.5*x(6)+1.9)/x(4))-1;
g(11)=((1.1*x(7)+1.9)/x(5))-1;

% No equality constraint in this problem, so empty;
geq=[];

% Apply inequality constraints
for k=1:length(g),
    Z=Z+ lam*g(k)^2*getH(g(k));
end
% Apply equality constraints
for k=1:length(geq),
   Z=Z+lam*geq(k)^2*getHeq(geq(k));
end
 obj=obj+Z;
end

% Four-bar truss problem
if strcmp(fun,'Bar')
    obj(1)=200*(2*x(1)+sqrt(2*x(2))+sqrt(x(3))+x(4));
    obj(2)=0.01*((2/x(2))+2*(sqrt(2)/x(2))-2*(sqrt(2)/x(3))+(2/x(4)));
end

% Disk brake design problem
if strcmp(fun,'Disk')
obj(1)=4.9*(10^(-5))*(x(2)^2-x(1)^2)*(x(4)-1);
obj(2)=9.82*10^6*(x(2)^2-x(1)^2)/((x(2)^3-x(1)^3)*x(4)*x(3));
Z=0;
% Penalty constant
lam=10^15;
g(1) = 20+x(1)-x(2);
g(2) = 2.5*(x(4)+1)-30;
g(3) = (x(3))/(3.14*(x(2)^2-x(1)^2)^2)-0.4;
g(4) = (2.22*10^(-3)*x(3)*(x(2)^3-x(1)^3))/((x(2)^2-x(1)^2)^2)-1;
g(5) = 900-(2.66*10^(-2)*x(3)*x(4)*(x(2)^3-x(1)^3))/((x(2)^2-x(1)^2));

% No equality constraint in this problem, so empty;
geq=[];

% Apply inequality constraints
for k=1:length(g),
    Z=Z+ lam*g(k)^2*getH(g(k));
end
% Apply equality constraints
for k=1:length(geq),
   Z=Z+lam*geq(k)^2*getHeq(geq(k));
end
obj=obj+Z;
end


end

function H=getH(g)
if g<=0, 
    H=0; 
else
    H=1; 
end
end
% Test if equalities hold
function H=geteqH(g)
if g==0,
    H=0;
else
    H=1; 
end
end


function Archive=DeleteFromRep(Archive,EXTRA,gamma)

    if nargin<3
        gamma=1;
    end

    for k=1:EXTRA
        [occ_cell_index occ_cell_member_count]=GetOccupiedCells(Archive);

        p=occ_cell_member_count.^gamma;
        p=p/sum(p);

        selected_cell_index=occ_cell_index(RouletteWheelSelection(p));

        GridIndices=[Archive.GridIndex];

        selected_cell_members=find(GridIndices==selected_cell_index);

        n=numel(selected_cell_members);

        selected_memebr_index=randi([1 n]);

        j=selected_cell_members(selected_memebr_index);
        
        Archive=[Archive(1:j-1); Archive(j+1:end)];
    end
    
end

3 运行结果

4 参考文献

[1] Ag A ,  Dz A ,  Am B . MOMRFO: Multi-objective Manta ray foraging optimizer for handling engineering design problems.  2021.

博主简介:擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真,相关matlab代码问题可私信交流。

部分理论引用网络文献,若有侵权联系博主删除。

  • 0
    点赞
  • 4
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

matlab科研助手

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值