From e46a5cdb9aae85124bee9c3d9e032e79e610945a Mon Sep 17 00:00:00 2001
From: Florentina2510 <97096352+Florentina2510@users.noreply.github.com>
Date: Tue, 4 Jan 2022 17:52:06 +0700
Subject: [PATCH] JOS (Joint Opposite Selection)

JOS: A well-match of Selective Leading Opposition (SLO) and Dynamic Opposite (DO)
---
 GWO_JOS.m | 254 +++++++++++++++++++++++++++++++++++++++++++++++++
 HHO_JOS.m | 280 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 MFO_JOS.m | 238 ++++++++++++++++++++++++++++++++++++++++++++++
 SOA_JOS.m | 213 +++++++++++++++++++++++++++++++++++++++++
 WOA_JOS.m | 244 +++++++++++++++++++++++++++++++++++++++++++++++
 5 files changed, 1229 insertions(+)
 create mode 100644 GWO_JOS.m
 create mode 100644 HHO_JOS.m
 create mode 100644 MFO_JOS.m
 create mode 100644 SOA_JOS.m
 create mode 100644 WOA_JOS.m

diff --git a/GWO_JOS.m b/GWO_JOS.m
new file mode 100644
index 0000000..ba5a6b9
--- /dev/null
+++ b/GWO_JOS.m
@@ -0,0 +1,254 @@
+function [Alpha_score,Alpha_pos,GWOJOS]=GWO_JOS(N,T,maxRun,maxFE,BFid,nD,fhd,Jr)
+
+if nargin ~= 8
+    N=30;               % Population size
+    maxRun = 10;        % Maximum Run
+    BFid = 1;           % Number id of benchmark function
+    nD = 10;            % Number of dimensions
+    maxFE = 10000*nD;   % Number of function evaluations
+    Jr=0.25;            % Jumping Rate
+    T=ceil(maxFE/N);    % Maximum number of iterations
+    fhd=str2func('cec17_func');
+end
+
+lb = -100*ones(1,nD); ub = 100*ones(1,nD);
+dim = nD;
+
+disp('GWO-JOS is now tackling your problem')
+
+GWOJOS=zeros(maxRun+1,T);
+
+for run=1:maxRun
+    tic
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %initializing boundary for opposition
+    Boundary_no= size(ub,2); % numnber of boundaries
+    
+    % If the boundaries of all variables are equal and user enter a signle
+    % number for both ub and lb
+    if Boundary_no==1
+        for i=1:dim
+            upper(1,i)=ub;
+            lower(1,i)=lb;
+        end
+        % If each variable has a different lb and ub
+    else
+        for i=1:dim
+            upper(1,i)=ub(i);
+            lower(1,i)=lb(i);
+        end
+    end
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
+    % initialize alpha, beta, and delta_pos
+    Alpha_pos=zeros(1,dim);
+    Alpha_score=inf; %change this to -inf for maximization problems
+    
+    Beta_pos=zeros(1,dim);
+    Beta_score=inf; %change this to -inf for maximization problems
+    
+    Delta_pos=zeros(1,dim);
+    Delta_score=inf; %change this to -inf for maximization problems
+    
+    %Initialize the positions of search agents
+    X=initialization(N,dim,ub,lb);
+    
+    for i=1:N
+        OP(i,:)=((ub-lb).*rand(size(lb)))+lb-X(i,:);
+        RO(i,:) = rand*OP(i,:);
+        DO(i,:) = X(i,:) + rand*(RO(i,:)-X(i,:));
+        
+        for j=1:dim
+            if DO(i,j)<lb(1,j)
+                DO(i,j)=lb(1,j);
+            end
+            if DO(i,j)>ub(1,j)
+                DO(i,j)=ub(1,j);
+            end
+        end
+    end
+    X=DO;
+    
+    t=0;% Loop counter
+    nFE=0;
+    Row_Id=1;
+    while nFE<maxFE
+        for i=1:size(X,1)
+            
+            % Return back the search agents that go beyond the boundaries of the search space
+            Flag4ub=X(i,:)>ub;
+            Flag4lb=X(i,:)<lb;
+            X(i,:)=(X(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;
+            
+            % Calculate objective function for each search agent
+            fitness(i)=feval(fhd,X(i,:)',BFid);
+            nFE = nFE + 1;
+            if nFE > maxFE; break; end
+            
+            % Update Alpha, Beta, and Delta
+            if fitness(i)<Alpha_score
+                Alpha_score=fitness(i); % Update alpha
+                Alpha_pos=X(i,:);
+                Row_Id=i;
+            end
+            
+            if fitness(i)>Alpha_score && fitness(i)<Beta_score
+                Beta_score=fitness(i); % Update beta
+                Beta_pos=X(i,:);
+                Row_Id=i;
+            end
+            
+            if fitness(i)>Alpha_score && fitness(i)>Beta_score && fitness(i)<Delta_score
+                Delta_score=fitness(i); % Update delta
+                Delta_pos=X(i,:);
+                Row_Id=i;
+            end
+            
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+            %updating boundary for opposition after every iteration
+            for x=1:size(X,1)
+                for y=1:size(X,2)
+                    if upper(1,y)<X(x,y)
+                        upper(1,y)=X(x,y);
+                    end
+                    if lower(1,y)>X(x,y)
+                        lower(1,y)=X(x,y);
+                    end
+                end
+            end
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        end
+        t=t+1;
+        GWOJOS(1,t) = nFE; %nFE
+        GWOJOS(run+1,t) = Alpha_score;
+        
+        a=2-t*((2)/T); % a decreases linearly fron 2 to 0
+        threshold=a;
+        X=corOppose2(X,ub,lb,upper,lower,dim,threshold,Row_Id);
+        
+        % Update the Position of search agents including omegas
+        for i=1:size(X,1)
+            for j=1:size(X,2)
+                
+                r1=rand(); % r1 is a random number in [0,1]
+                r2=rand(); % r2 is a random number in [0,1]
+                
+                A1=2*a*r1-a; % Equation (3.3)
+                C1=2*r2; % Equation (3.4)
+                
+                D_alpha=abs(C1*Alpha_pos(j)-X(i,j)); % Equation (3.5)-part 1
+                X1=Alpha_pos(j)-A1*D_alpha; % Equation (3.6)-part 1
+                
+                r1=rand();
+                r2=rand();
+                
+                A2=2*a*r1-a; % Equation (3.3)
+                C2=2*r2; % Equation (3.4)
+                
+                D_beta=abs(C2*Beta_pos(j)-X(i,j)); % Equation (3.5)-part 2
+                X2=Beta_pos(j)-A2*D_beta; % Equation (3.6)-part 2
+                
+                r1=rand();
+                r2=rand();
+                
+                A3=2*a*r1-a; % Equation (3.3)
+                C3=2*r2; % Equation (3.4)
+                
+                D_delta=abs(C3*Delta_pos(j)-X(i,j)); % Equation (3.5)-part 3
+                X3=Delta_pos(j)-A3*D_delta; % Equation (3.5)-part 3
+                
+                X(i,j)=(X1+X2+X3)/3;% Equation (3.7)
+                
+            end
+        end
+        
+        if rand < Jr && nFE+N < maxFE
+            
+            for i=1:N
+                OP(i,:)=((ub-lb).*rand(size(lb)))+lb-X(i,:);
+                RO(i,:) = rand*OP(i,:);
+                DO(i,:) = X(i,:) + rand*(RO(i,:)-X(i,:));
+                for j=1:dim
+                    if DO(i,j)<lb(1,j)
+                        DO(i,j)=lb(1,j);
+                    end
+                    if DO(i,j)>ub(1,j)
+                        DO(i,j)=ub(1,j);
+                    end
+                end
+            end
+            X=DO;
+        end
+    end
+    toc
+end
+% display(['The best solution obtained by GWO is : ', num2str(Best_pos)]);
+% display(['The best optimal value of the objective funciton found by GWO is : ', num2str(Best_score)]);
+end
+
+function Positions=initialization(N,dim,ub,lb)
+Boundary_no= size(ub,2); % numnber of boundaries
+
+% If the boundaries of all variables are equal and user enter a signle
+% number for both ub and lb
+if Boundary_no==1
+    Positions=rand(N,dim).*(ub-lb)+lb;
+end
+
+% If each variable has a different lb and ub
+if Boundary_no>1
+    for i=1:dim
+        ub_i=ub(i);
+        lb_i=lb(i);
+        Positions(:,i)=rand(N,1).*(ub_i-lb_i)+lb_i;
+    end
+end
+end
+
+function [Positions]=corOppose2(Positions,ub,lb,upper,lower,dim,threshold,Row_Id)
+
+[n,b] = size(Positions);
+for i=1:n(1) 
+    
+    if i ~= Row_Id
+        
+        sum=0;
+        greater=[];
+        less=[];
+        x=1;z=1;y=1;
+        
+        for j=1:b
+            d(x)=abs(Positions(Row_Id,j)-Positions(i,j));
+            if d(x)<threshold
+                greater(y)=j;
+                y=y+1;
+            else
+                less(z)=j;
+                z=z+1;
+            end
+            sum=sum+d(x)*d(x);
+            x=x+1;
+        end
+        %     sum
+        src=1-(double(6*sum))/(double(n(1)*(n(1)*n(1)-1)));
+        %     src
+        if src<=0
+            if size(greater)<size(less)
+                %             for j=1:size(less)
+                %                 dim=less(j);
+                %                 Positions(i,dim)=ub(dim)+lb(dim)-Positions(i,dim);
+                %             end
+            else
+                for j=1:size(greater)
+                    dim=greater(j);
+                    Positions(i,dim)=(upper(1,dim)+lower(1,dim)-Positions(i,dim));
+                end
+            end
+        end
+    end
+end
+end
+
+
+
+
diff --git a/HHO_JOS.m b/HHO_JOS.m
new file mode 100644
index 0000000..8b17fe3
--- /dev/null
+++ b/HHO_JOS.m
@@ -0,0 +1,280 @@
+function [Rabbit_Energy,Rabbit_Location,HHOJOS]=HHO_JOS(N,T,maxRun,maxFE,BFid,nD,fhd,Jr)
+
+if nargin ~= 8
+    N=30;               % Population size
+    maxRun = 10;        % Maximum Run
+    BFid = 1;           % Number id of benchmark function
+    nD = 10;            % Number of dimensions
+    maxFE = 10000*nD;   % Number of function evaluations
+    Jr=0.25;            % Jumping Rate
+    T=ceil(maxFE/N);    % Maximum number of iterations
+    fhd=str2func('cec17_func');
+end
+
+lb = -100*ones(1,nD); ub = 100*ones(1,nD);
+dim = nD;
+
+disp('HHO-JOS is now tackling your problem')
+
+HHOJOS=zeros(maxRun+1,T);
+
+for run=1:maxRun
+    tic
+    % initialize the location and Energy of the rabbit
+    Rabbit_Location=zeros(1,dim);
+    Rabbit_Energy=inf;
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %initializing boundary for opposition
+    Boundary_no= size(ub,2); % numnber of boundaries
+    
+    % If the boundaries of all variables are equal and user enter a signle
+    % number for both ub and lb
+    if Boundary_no==1
+        for i=1:dim
+            upper(1,i)=ub;
+            lower(1,i)=lb;
+        end
+        % If each variable has a different lb and ub
+    else
+        for i=1:dim
+            upper(1,i)=ub(i);
+            lower(1,i)=lb(i);
+        end
+    end
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
+    %Initialize the locations of Harris' hawks
+    X=initialization(N,dim,ub,lb);
+    %X=initialization(N,dim,up,down);
+    for i=1:N
+        OP(i,:)=((ub-lb).*rand(size(lb)))+lb-X(i,:);
+        RO(i,:) = rand*OP(i,:);
+        DO(i,:) = X(i,:) + rand*(RO(i,:)-X(i,:));
+        
+        for j=1:dim
+            if DO(i,j)<lb(1,j)
+                DO(i,j)=lb(1,j);
+            end
+            if DO(i,j)>ub(1,j)
+                DO(i,j)=ub(1,j);
+            end
+        end
+    end
+    X=DO;
+    
+    t=0; % Loop counter
+    nFE=0;
+    h=size(X,1);
+    fitness = zeros(1,h); %row
+
+    while nFE<maxFE
+ 
+        for i=1:h
+            % Check boundries
+            FU=X(i,:)>ub;
+            FL=X(i,:)<lb;
+            X(i,:)=(X(i,:).*(~(FU+FL)))+ub.*FU+lb.*FL;
+
+            fitness(i)=feval(fhd,X(i,:)',BFid);
+            nFE = nFE + 1;
+            if nFE > maxFE; break; end        
+            
+            % Update the location of Rabbit
+            if fitness(i)<Rabbit_Energy
+                Rabbit_Energy=fitness(i);
+                Rabbit_Location=X(i,:);
+                Rabbit_Row_Id = i;  %%
+            end
+            
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+            %updating boundary for opposition after every iteration
+            for x=1:size(X,1)
+                for y=1:size(X,2)
+                    if upper(1,y)<X(x,y)
+                        upper(1,y)=X(x,y);
+                    end
+                    if lower(1,y)>X(x,y)
+                        lower(1,y)=X(x,y);
+                    end
+                end
+            end
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        end
+        
+        %nFE=nFE+h;
+        t=t+1;
+        HHOJOS(1,t) = nFE; %nFE
+        HHOJOS(run+1,t) = Rabbit_Energy;
+        
+        E1=2*(1-(t/T)); % factor to show the decreaing energy of rabbit
+        
+        % Oppose the least fitness elements
+        threshold=E1;
+      
+        X=corOppose2(X,fitness,ub,lb,upper,lower,dim,threshold,Rabbit_Row_Id);
+        
+        % Update the location of Harris' hawks
+        for i=1:h
+            E0=2*rand()-1; %-1<E0<1
+            Escaping_Energy=E1*(E0);  % escaping energy of rabbit
+            
+            if abs(Escaping_Energy)>=1 %15.4 percent
+                %% Exploration:
+                % Harris' hawks perch randomly based on 2 strategy:
+                
+                q=rand();
+                rand_Hawk_index = floor(N*rand()+1);
+                X_rand = X(rand_Hawk_index, :);
+                if q<0.5
+                    % perch based on other family members
+                    X(i,:)=X_rand-rand()*abs(X_rand-2*rand()*X(i,:)); %TLBO
+                elseif q>=0.5
+                    % perch on a random tall tree (random site inside group's home range)
+                    X(i,:)=(Rabbit_Location(1,:)-mean(X))-rand()*((ub-lb)*rand+lb); %QPSO
+                end
+                
+            elseif abs(Escaping_Energy)<1
+                %% Exploitation: 84.6 percent
+                % 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  %%21.149 percent
+                    X(i,:)=(Rabbit_Location)-Escaping_Energy*abs(Rabbit_Location-X(i,:));
+                end
+                
+                if r>=0.5 && abs(Escaping_Energy)>=0.5  % Soft besiege %% 21.149
+                    Jump_strength=2*(1-rand()); % random jump strength of the rabbit
+                    X(i,:)=(Rabbit_Location-X(i,:))-Escaping_Energy*abs(Jump_strength*Rabbit_Location-X(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()); % 0< Jump_strength <2
+                    X1=Rabbit_Location-Escaping_Energy*abs(Jump_strength*Rabbit_Location-X(i,:));
+                    
+                    if feval(fhd,X1',BFid) < feval(fhd,X(i,:)',BFid)
+                        X(i,:)=X1;
+                        
+                    else % hawks perform levy-based short rapid dives around the rabbit
+                        X2=Rabbit_Location-Escaping_Energy*abs(Jump_strength*Rabbit_Location-X(i,:))+rand(1,dim).*Levy(dim);
+                        if feval(fhd,X2',BFid) < feval(fhd,X(i,:)',BFid)
+                            X(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(X));
+                    
+                    if feval(fhd,X1',BFid) < feval(fhd,X(i,:)',BFid)
+                        X(i,:)=X1;
+                    else % Perform levy-based short rapid dives around the rabbit
+                        X2=Rabbit_Location-Escaping_Energy*abs(Jump_strength*Rabbit_Location-mean(X))+rand(1,dim).*Levy(dim);
+                        if feval(fhd,X2',BFid) < feval(fhd,X(i,:)',BFid)
+                            X(i,:)=X2;
+                        end
+                    end
+                end
+                %%
+            end
+        end
+        
+        
+        if rand < Jr && nFE+h < maxFE
+            
+            for i=1:N
+                OP(i,:)=((ub-lb).*rand(size(lb)))+lb-X(i,:);
+                RO(i,:) = rand*OP(i,:);
+                DO(i,:) = X(i,:) + rand*(RO(i,:)-X(i,:));
+                for j=1:dim
+                    if DO(i,j)<lb(1,j)
+                        DO(i,j)=lb(1,j);
+                    end
+                    if DO(i,j)>ub(1,j)
+                        DO(i,j)=ub(1,j);
+                    end
+                end
+            end
+            X=DO;
+        end
+    end
+    toc
+end
+end
+
+% ___________________________________
+function o=Levy(d)
+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,d)*sigma;v=randn(1,d);
+step=u./abs(v).^(1/beta);
+o=step;
+end
+
+function [X]=initialization(N,dim,up,down)
+
+if size(up,1)==1
+    X=rand(N,dim).*repmat((up-down)+down,N,1);
+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 [Positions]=corOppose2(Positions,fitness,ub,lb,upper,lower,dim,threshold,Rabbit_Row_Id)
+
+[n,b] = size(Positions);
+for i=1:n(1) 
+    
+    if i ~= Rabbit_Row_Id
+        
+        sum=0;
+        greater=[];
+        less=[];
+        x=1;z=1;y=1;
+        
+        for j=1:b
+            d(x)=abs(Positions(Rabbit_Row_Id,j)-Positions(i,j));
+            if d(x)<threshold
+                greater(y)=j;
+                y=y+1;
+            else
+                less(z)=j;
+                z=z+1;
+            end
+            sum=sum+d(x)*d(x);
+            x=x+1;
+        end
+        %     sum
+        src=1-(double(6*sum))/(double(n(1)*(n(1)*n(1)-1)));
+        %     src
+        if src<=0
+            if size(greater)<size(less)
+                %             for j=1:size(less)
+                %                 dim=less(j);
+                %                 Positions(i,dim)=ub(dim)+lb(dim)-Positions(i,dim);
+                %             end
+            else
+                for j=1:size(greater)
+                    dim=greater(j);
+                    Positions(i,dim)=(upper(1,dim)+lower(1,dim)-Positions(i,dim));
+                end
+            end
+        end
+    end
+end
+end
+
+
+
diff --git a/MFO_JOS.m b/MFO_JOS.m
new file mode 100644
index 0000000..fd4d24d
--- /dev/null
+++ b/MFO_JOS.m
@@ -0,0 +1,238 @@
+function [Best_flame_score,Best_flame_pos,MFOJOS]=MFO_JOS(N,T,maxRun,maxFE,BFid,nD,fhd,Jr)
+
+if nargin ~= 8
+    N=30;               % Population size
+    maxRun = 10;        % Maximum Run
+    BFid = 1;           % Number id of benchmark function
+    nD = 10;            % Number of dimensions
+    maxFE = 10000*nD;   % Number of function evaluations
+    Jr=0.25;            % Jumping Rate
+    T=ceil(maxFE/N);    % Maximum number of iterations
+    fhd=str2func('cec17_func');
+end
+
+lb=-100*ones(1,nD); ub=100*ones(1,nD);
+dim = nD;
+
+MFOJOS=zeros(maxRun+1,T);
+
+for run=1:maxRun
+    tic
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %initializing boundary for opposition
+    Boundary_no= size(ub,2); % numnber of boundaries
+    
+    % If the boundaries of all variables are equal and user enter a signle
+    % number for both ub and lb
+    if Boundary_no==1
+        for i=1:dim
+            upper(1,i)=ub;
+            lower(1,i)=lb;
+        end
+        % If each variable has a different lb and ub
+    else
+        for i=1:dim
+            upper(1,i)=ub(i);
+            lower(1,i)=lb(i);
+        end
+    end
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %Initialize the positions of moths
+    Moth_pos=initialization(N,dim,ub,lb);
+    
+    for i=1:N
+        OP(i,:)=((ub-lb).*rand(size(lb)))+lb-Moth_pos(i,:);
+        RO(i,:) = rand*OP(i,:);
+        DO(i,:) = Moth_pos(i,:) + rand*(RO(i,:)-Moth_pos(i,:));
+        
+        for j=1:dim
+            if DO(i,j)<lb(1,j)
+                DO(i,j)=lb(1,j);
+            end
+            if DO(i,j)>ub(1,j)
+                DO(i,j)=ub(1,j);
+            end
+        end
+    end
+    Moth_pos=DO;
+    
+    nFE=0;
+    Iter=0;
+    Row_Id=1;
+    while nFE<maxFE
+        % Number of flames Eq. (3.14) in the paper MFO
+        Flame_no=round(N-Iter*((N-1)/T));
+        for i=1:size(Moth_pos,1)
+            % Check if moths go out of the search spaceand bring it back
+            Flag4ub=Moth_pos(i,:)>ub;
+            Flag4lb=Moth_pos(i,:)<lb;
+            Moth_pos(i,:)=(Moth_pos(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;
+            % Calculate the fitness of moths
+            Moth_fitness(1,i)=feval(fhd,Moth_pos(i,:)',BFid);
+            nFE=nFE+1;
+         
+            if Iter==0
+                % Sort the first population of moths
+                [fitness_sorted I]=sort(Moth_fitness);
+                sorted_population=Moth_pos(I,:);
+                % Update the flames
+                best_flames=sorted_population;
+                best_flame_fitness=fitness_sorted;
+                Row_Id=i;
+            else
+                % Sort the moths
+                double_population=[previous_population;best_flames];
+                double_fitness=[previous_fitness best_flame_fitness];
+                
+                [double_fitness_sorted, I]=sort(double_fitness);
+                double_sorted_population=double_population(I,:);
+                
+                fitness_sorted=double_fitness_sorted(1:N);
+                sorted_population=double_sorted_population(1:N,:);
+                
+                % Update the flames
+                best_flames=sorted_population;
+                best_flame_fitness=fitness_sorted;
+                Row_Id=i;
+            end
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+            %updating boundary for opposition after every iteration
+            for x=1:size(Moth_pos,1)
+                for y=1:size(Moth_pos,2)
+                    if upper(1,y)<Moth_pos(x,y)
+                        upper(1,y)=Moth_pos(x,y);
+                    end
+                    if lower(1,y)>Moth_pos(x,y)
+                        lower(1,y)=Moth_pos(x,y);
+                    end
+                end
+            end
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        end
+        
+        % Update the position best flame obtained so far
+        Best_flame_score=fitness_sorted(1);
+        Best_flame_pos=sorted_population(1,:);
+        
+        previous_population=Moth_pos;
+        previous_fitness=Moth_fitness;
+        
+        Iter=Iter+1;
+        MFOJOS(1,Iter)=nFE;
+        MFOJOS(run+1,Iter)=Best_flame_score;
+        
+        % a linearly dicreases from -1 to -2 to calculate t in Eq. (3.12)
+        a=-1+Iter*((-1)/T);
+        threshold=a;
+        
+        Moth_pos=corOppose2(Moth_pos,ub,lb,upper,lower,dim,threshold,Row_Id);
+        
+        for i=1:size(Moth_pos,1)
+            for j=1:size(Moth_pos,2)
+                if i<=Flame_no % Update the position of the moth with respect to its corresponsing flame
+                    % D in Eq. (3.13) in MFO paper
+                    distance_to_flame=abs(sorted_population(i,j)-Moth_pos(i,j));
+                    b=1;
+                    t=(a-1)*rand+1;
+                    % Eq. (3.12) in MFO paper
+                    Moth_pos(i,j)=distance_to_flame*exp(b.*t).*cos(t.*2*pi)+sorted_population(i,j);
+                end
+                
+                if i>Flame_no % Upaate the position of the moth with respct to one flame
+                    % Eq. (3.13) in MFO paper
+                    distance_to_flame=abs(sorted_population(i,j)-Moth_pos(i,j));
+                    b=1;
+                    t=(a-1)*rand+1;
+                    % Eq. (3.12) in MFO paper
+                    Moth_pos(i,j)=distance_to_flame*exp(b.*t).*cos(t.*2*pi)+sorted_population(Flame_no,j);
+                end
+            end
+        end
+        
+        if rand < Jr && nFE+N < maxFE
+            
+            for i=1:N
+                OP(i,:)=((ub-lb).*rand(size(lb)))+lb-Moth_pos(i,:);
+                RO(i,:) = rand*OP(i,:);
+                DO(i,:) = Moth_pos(i,:) + rand*(RO(i,:)-Moth_pos(i,:));
+                for j=1:dim
+                    if DO(i,j)<lb(1,j)
+                        DO(i,j)=lb(1,j);
+                    end
+                    if DO(i,j)>ub(1,j)
+                        DO(i,j)=ub(1,j);
+                    end
+                end
+            end
+            Moth_pos=DO;
+        end
+    end
+    toc
+end
+% display(['The best solution obtained by MFO is : ', num2str(Best_flame_score)]);
+% display(['The best optimal value of the objective funciton found by MFO is : ', num2str(Best_flame_pos)]);
+end
+
+function X=initialization(SearchAgents_no,dim,ub,lb)
+Boundary_no= size(ub,2); % numnber of boundaries
+% If the boundaries of all variables are equal and user enter a signle
+% number for both ub and lb
+if Boundary_no==1
+    X=rand(SearchAgents_no,dim).*(ub-lb)+lb;
+end
+% If each variable has a different lb and ub
+if Boundary_no>1
+    for i=1:dim
+        ub_i=ub(i);
+        lb_i=lb(i);
+        X(:,i)=rand(SearchAgents_no,1).*(ub_i-lb_i)+lb_i;
+    end
+end
+end
+
+function [Positions]=corOppose2(Positions,ub,lb,upper,lower,dim,threshold,Row_Id)
+
+[n,b] = size(Positions);
+for i=1:n(1) 
+    
+    if i ~= Row_Id
+        
+        sum=0;
+        greater=[];
+        less=[];
+        x=1;z=1;y=1;
+        
+        for j=1:b
+            d(x)=abs(Positions(Row_Id,j)-Positions(i,j));
+            if d(x)<threshold
+                greater(y)=j;
+                y=y+1;
+            else
+                less(z)=j;
+                z=z+1;
+            end
+            sum=sum+d(x)*d(x);
+            x=x+1;
+        end
+        %     sum
+        src=1-(double(6*sum))/(double(n(1)*(n(1)*n(1)-1)));
+        %     src
+        if src<=0
+            if size(greater)<size(less)
+                %             for j=1:size(less)
+                %                 dim=less(j);
+                %                 Positions(i,dim)=ub(dim)+lb(dim)-Positions(i,dim);
+                %             end
+            else
+                for j=1:size(greater)
+                    dim=greater(j);
+                    Positions(i,dim)=(upper(1,dim)+lower(1,dim)-Positions(i,dim));
+                end
+            end
+        end
+    end
+end
+end
+
+
+
diff --git a/SOA_JOS.m b/SOA_JOS.m
new file mode 100644
index 0000000..b743075
--- /dev/null
+++ b/SOA_JOS.m
@@ -0,0 +1,213 @@
+function [Score,Position,SOAJOS]=SOA_JOS(N,T,maxRun,maxFE,BFid,nD,fhd,Jr)
+
+if nargin ~= 8
+    N=30;               % Population size
+    maxRun = 10;        % Maximum Run
+    BFid = 1;           % Number id of benchmark function
+    nD = 10;            % Number of dimensions
+    maxFE = 10000*nD;   % Number of function evaluations
+    Jr=0.25;            % Jumping Rate
+    T=ceil(maxFE/N);    % Maximum number of iterations
+    fhd=str2func('cec17_func');
+end
+
+lb=-100*ones(1,nD); ub=100*ones(1,nD);
+dim = nD;
+
+disp('SOA-JOS is now tackling your problem')
+
+SOAJOS=zeros(maxRun+1,T);
+
+for run=1:maxRun
+    tic
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %initializing boundary for opposition
+    Boundary_no= size(ub,2); % numnber of boundaries
+    
+    % If the boundaries of all variables are equal and user enter a signle
+    % number for both ub and lb
+    if Boundary_no==1
+        for i=1:dim
+            upper(1,i)=ub;
+            lower(1,i)=lb;
+        end
+        % If each variable has a different lb and ub
+    else
+        for i=1:dim
+            upper(1,i)=ub(i);
+            lower(1,i)=lb(i);
+        end
+    end
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
+    
+    Position=zeros(1,dim);
+    Score=inf;
+    
+    Positions=init(N,dim,ub,lb);
+    
+    for i=1:N
+        OP(i,:)=((ub-lb).*rand(size(lb)))+lb-Positions(i,:);
+        RO(i,:) = rand*OP(i,:);
+        DO(i,:) = Positions(i,:) + rand*(RO(i,:)-Positions(i,:));
+        
+        for j=1:dim
+            if DO(i,j)<lb(1,j)
+                DO(i,j)=lb(1,j);
+            end
+            if DO(i,j)>ub(1,j)
+                DO(i,j)=ub(1,j);
+            end
+        end
+    end
+    Positions=DO;
+    
+    t=0;
+    nFE=0;
+    Row_Id=1;
+    
+    while nFE<maxFE
+        
+        for i=1:size(Positions,1)
+            
+            Flag4ub=Positions(i,:)>ub;
+            Flag4lb=Positions(i,:)<lb;
+            Positions(i,:)=(Positions(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;
+            
+            fitness=feval(fhd,Positions(i,:)',BFid);
+            nFE = nFE + 1;
+            if nFE > maxFE; break; end
+            if fitness<Score
+                Score=fitness;
+                Position=Positions(i,:);
+                Row_Id=i;
+            end
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+            %updating boundary for opposition after every iteration
+            for x=1:size(Positions,1)
+                for y=1:size(Positions,2)
+                    if upper(1,y)<Positions(x,y)
+                        upper(1,y)=Positions(x,y);
+                    end
+                    if lower(1,y)>Positions(x,y)
+                        lower(1,y)=Positions(x,y);
+                    end
+                end
+            end
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        end
+        
+        t=t+1;
+        SOAJOS(1,t) = nFE; %nFE
+        SOAJOS(run+1,t) = Score;
+        
+        Fc=2-t*((2)/T);
+        
+        % Oppose the least fitness elements
+        threshold=Fc;
+        
+        Positions=corOppose2(Positions,ub,lb,upper,lower,dim,threshold,Row_Id);
+        
+        for i=1:size(Positions,1)
+            for j=1:size(Positions,2)
+                
+                r1=rand();
+                r2=rand();
+                
+                A1=2*Fc*r1-Fc;
+                C1=2*r2;
+                b=1;
+                ll=(Fc-1)*rand()+1;
+                
+                D_alphs=Fc*Positions(i,j)+A1*((Position(j)-Positions(i,j)));
+                X1=D_alphs*exp(b.*ll).*cos(ll.*2*pi)+Position(j);
+                Positions(i,j)=X1;
+                
+            end
+        end
+        if rand < Jr && nFE+size(Positions,1) < maxFE
+            
+            for i=1:N
+                OP(i,:)=((ub-lb).*rand(size(lb)))+lb-Positions(i,:);
+                RO(i,:) = rand*OP(i,:);
+                DO(i,:) = Positions(i,:) + rand*(RO(i,:)-Positions(i,:));
+                for j=1:dim
+                    if DO(i,j)<lb(1,j)
+                        DO(i,j)=lb(1,j);
+                    end
+                    if DO(i,j)>ub(1,j)
+                        DO(i,j)=ub(1,j);
+                    end
+                end
+            end
+            Positions=DO;
+        end
+    end
+    % display(['The best solution obtained by SOA is : ', num2str(Position)]);
+    % display(['The best optimal value of the objective funciton found by SOA is : ', num2str(Score)]);
+    toc
+end
+end
+function Pos=init(SearchAgents,dim,ub,lb)
+
+Boundary= size(ub,2);
+if Boundary==1
+    Pos=rand(SearchAgents,dim).*(ub-lb)+lb;
+end
+
+if Boundary>1
+    for i=1:dim
+        ub_i=ub(i);
+        lb_i=lb(i);
+        Pos(:,i)=rand(SearchAgents,1).*(ub_i-lb_i)+lb_i;
+    end
+end
+end
+
+function [Positions]=corOppose2(Positions,ub,lb,upper,lower,dim,threshold,Row_Id)
+
+[n,b] = size(Positions);
+for i=1:n(1) 
+    
+    if i ~= Row_Id
+        
+        sum=0;
+        greater=[];
+        less=[];
+        x=1;z=1;y=1;
+        
+        for j=1:b
+            d(x)=abs(Positions(Row_Id,j)-Positions(i,j));
+            if d(x)<threshold
+                greater(y)=j;
+                y=y+1;
+            else
+                less(z)=j;
+                z=z+1;
+            end
+            sum=sum+d(x)*d(x);
+            x=x+1;
+        end
+        %     sum
+        src=1-(double(6*sum))/(double(n(1)*(n(1)*n(1)-1)));
+        %     src
+        if src<=0
+            if size(greater)<size(less)
+                %             for j=1:size(less)
+                %                 dim=less(j);
+                %                 Positions(i,dim)=ub(dim)+lb(dim)-Positions(i,dim);
+                %             end
+            else
+                for j=1:size(greater)
+                    dim=greater(j);
+                    Positions(i,dim)=(upper(1,dim)+lower(1,dim)-Positions(i,dim));
+                end
+            end
+        end
+    end
+end
+end
+
+
+
diff --git a/WOA_JOS.m b/WOA_JOS.m
new file mode 100644
index 0000000..319644f
--- /dev/null
+++ b/WOA_JOS.m
@@ -0,0 +1,244 @@
+function [Leader_score,Leader_pos,WOACNVG]=WOA_JOS(N,T,maxRun,maxFE,BFid,nD,fhd,Jr)
+
+if nargin ~= 8
+    N=30;               % Population size
+    maxRun = 10;        % Maximum Run
+    BFid = 1;           % Number id of benchmark function
+    nD = 10;            % Number of dimensions
+    maxFE = 10000*nD;   % Number of function evaluations
+    Jr=0.25;            % Jumping Rate
+    T=ceil(maxFE/N);    % Maximum number of iterations
+    fhd=str2func('cec17_func');
+end
+
+lb=-100*ones(1,nD); ub=100*ones(1,nD);
+dim = nD;
+
+disp('WOA-JOS is now tackling your problem')
+
+WOACNVG=zeros(maxRun+1,T);
+
+for run=1:maxRun
+    tic
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %initializing boundary for opposition
+    Boundary_no= size(ub,2); % numnber of boundaries
+    
+    % If the boundaries of all variables are equal and user enter a signle
+    % number for both ub and lb
+    if Boundary_no==1
+        for i=1:dim
+            upper(1,i)=ub;
+            lower(1,i)=lb;
+        end
+        % If each variable has a different lb and ub
+    else
+        for i=1:dim
+            upper(1,i)=ub(i);
+            lower(1,i)=lb(i);
+        end
+    end
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
+    % initialize position vector and score for the leader
+    Leader_pos=zeros(1,dim);
+    Leader_score=inf; %change this to -inf for maximization problems
+    
+    %Initialize the positions of search agents
+    X=initialization(N,dim,ub,lb);
+    
+    for i=1:N
+        OP(i,:)=((ub-lb).*rand(size(lb)))+lb-X(i,:);
+        RO(i,:) = rand*OP(i,:);
+        DO(i,:) = X(i,:) + rand*(RO(i,:)-X(i,:));
+        
+        for j=1:dim
+            if DO(i,j)<lb(1,j)
+                DO(i,j)=lb(1,j);
+            end
+            if DO(i,j)>ub(1,j)
+                DO(i,j)=ub(1,j);
+            end
+        end
+    end
+    X=DO;
+    
+    t=0;% Loop counter
+    nFE=0;
+    Row_Id=1;
+    
+    while nFE<maxFE
+        for i=1:size(X,1)
+            
+            % Return back the search agents that go beyond the boundaries of the search space
+            Flag4ub=X(i,:)>ub;
+            Flag4lb=X(i,:)<lb;
+            X(i,:)=(X(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;
+            
+            % Calculate objective function for each search agent
+            %fitness=fobj(X(i,:));
+            
+            fitness=feval(fhd,X(i,:)',BFid);
+            nFE = nFE + 1;
+            if nFE > maxFE; break; end
+            
+            % Update the leader
+            if fitness<Leader_score % Change this to > for maximization problem
+                Leader_score=fitness; % Update alpha
+                Leader_pos=X(i,:);
+                Row_Id = i;
+            end
+            
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+            %updating boundary for opposition after every iteration
+            for x=1:size(X,1)
+                for y=1:size(X,2)
+                    if upper(1,y)<X(x,y)
+                        upper(1,y)=X(x,y);
+                    end
+                    if lower(1,y)>X(x,y)
+                        lower(1,y)=X(x,y);
+                    end
+                end
+            end
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        end
+        t=t+1;
+        WOACNVG(1,t) = nFE; %nFE
+        WOACNVG(run+1,t) = Leader_score;
+        
+        a=2-t*((2)/T); % a decreases linearly fron 2 to 0
+        % Oppose the least fitness elements
+        threshold=a;
+        
+        X=corOppose2(X,ub,lb,upper,lower,dim,threshold,Row_Id);
+        % a2 linearly dicreases from -1 to -2 to calculate t in Eq. (3.12) in WOA paper
+
+        a2=-1+t*((-1)/T);
+        
+        % Update the Position of search agents
+        for i=1:size(X,1)
+            r1=rand(); % r1 is a random number in [0,1] in WOA paper
+            r2=rand(); % r2 is a random number in [0,1] in WOA paper
+            
+            A=2*a*r1-a;  % Eq. (2.3) in the paper of WOA
+            C=2*r2;      % Eq. (2.4) in the paper of WOA
+            
+            
+            b=1;               %  parameters in Eq. (2.5) in WOA paper
+            l=(a2-1)*rand+1;   %  parameters in Eq. (2.5) in WOA paper
+            
+            p = rand();        % p in Eq. (2.6) in WOA paper
+            
+            for j=1:size(X,2)
+                
+                if p<0.5
+                    if abs(A)>=1
+                        rand_leader_index = floor(N*rand()+1);
+                        X_rand = X(rand_leader_index, :);
+                        D_X_rand=abs(C*X_rand(j)-X(i,j)); % Eq. (2.7) in WOA paper
+                        X(i,j)=X_rand(j)-A*D_X_rand;      % Eq. (2.8) in WOA paper
+                        
+                    elseif abs(A)<1
+                        D_Leader=abs(C*Leader_pos(j)-X(i,j)); % Eq. (2.1) in WOA paper
+                        X(i,j)=Leader_pos(j)-A*D_Leader;      % Eq. (2.2) in WOA paper
+                    end
+                    
+                elseif p>=0.5
+                    
+                    distance2Leader=abs(Leader_pos(j)-X(i,j));
+                    % Eq. (2.5) in WOA paper
+                    X(i,j)=distance2Leader*exp(b.*l).*cos(l.*2*pi)+Leader_pos(j);
+                end
+            end
+        end
+        
+        if rand < Jr && nFE+size(X,1) < maxFE
+            
+            for i=1:N
+                OP(i,:)=((ub-lb).*rand(size(lb)))+lb-X(i,:);
+                RO(i,:) = rand*OP(i,:);
+                DO(i,:) = X(i,:) + rand*(RO(i,:)-X(i,:));
+                for j=1:dim
+                    if DO(i,j)<lb(1,j)
+                        DO(i,j)=lb(1,j);
+                    end
+                    if DO(i,j)>ub(1,j)
+                        DO(i,j)=ub(1,j);
+                    end
+                end
+            end
+            X=DO;
+        end
+    end
+    toc
+end
+% display(['The best solution obtained by GWO is : ', num2str(Best_pos)]);
+% display(['The best optimal value of the objective funciton found by GWO is : ', num2str(Best_score)]);
+end
+
+function Positions=initialization(N,dim,ub,lb)
+Boundary_no= size(ub,2); % numnber of boundaries
+
+% If the boundaries of all variables are equal and user enter a signle
+% number for both ub and lb
+if Boundary_no==1
+    Positions=rand(N,dim).*(ub-lb)+lb;
+end
+
+% If each variable has a different lb and ub
+if Boundary_no>1
+    for i=1:dim
+        ub_i=ub(i);
+        lb_i=lb(i);
+        Positions(:,i)=rand(N,1).*(ub_i-lb_i)+lb_i;
+    end
+end
+end
+
+function [Positions]=corOppose2(Positions,ub,lb,upper,lower,dim,threshold,Row_Id)
+
+[n,b] = size(Positions);
+for i=1:n(1) 
+    
+    if i ~= Row_Id
+        
+        sum=0;
+        greater=[];
+        less=[];
+        x=1;z=1;y=1;
+        
+        for j=1:b
+            d(x)=abs(Positions(Row_Id,j)-Positions(i,j));
+            if d(x)<threshold
+                greater(y)=j;
+                y=y+1;
+            else
+                less(z)=j;
+                z=z+1;
+            end
+            sum=sum+d(x)*d(x);
+            x=x+1;
+        end
+        %     sum
+        src=1-(double(6*sum))/(double(n(1)*(n(1)*n(1)-1)));
+        %     src
+        if src<=0
+            if size(greater)<size(less)
+                %             for j=1:size(less)
+                %                 dim=less(j);
+                %                 Positions(i,dim)=ub(dim)+lb(dim)-Positions(i,dim);
+                %             end
+            else
+                for j=1:size(greater)
+                    dim=greater(j);
+                    Positions(i,dim)=(upper(1,dim)+lower(1,dim)-Positions(i,dim));
+                end
+            end
+        end
+    end
+end
+end
+
+
+