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

原创已于 2022-08-02 16:22:01 修改 · 409 阅读

5 ·

CC 4.0 BY-SA版权

文章标签：

#算法 #matlab #开发语言

于 2022-07-30 09:25:08 首次发布

优化求解专栏收录该内容

553 篇文章

订阅专栏

该文研究了一种新的多目标优化算法——多目标蝠鲼觅食优化器（MOMRFO），用于解决工程设计中的多目标问题。MOMRFO结合了精英概念和外部存档，以维护种群的多样性和收敛性。通过仿真多个测试函数和实际工程问题，实验结果表明MOMRFO在收敛性和多样性上表现出色，证明了其在处理现实世界多目标问题上的潜力。

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