function [Ceq1_fit,valuess, Ceq1]=MS_EO(fhd,LB,UB,D,iter,swarm_size,varargin)
%Copyright (c) 2021, All rights reserved.
%Everyone is permitted to copy and distribute verbatim copies
%of this license document, but changing it is not allowed.

%  Mulit-strategy Synthetized Equilibrium Optimizer, MS-EO
%  Developed in MATLAB R2014b                                             %
%  Programmer: Qiuying Ling and Xinming Zhang                             %
%  E-Mail: xinmingzhang@126.com                                           %
%  Date: 2022-5-23
% Input:
% fhd; objective function
% LB: the lower boundary vector of the search space
% UB: the upper boundary vector of the search sapce
% D: the dimension of the search sapce
% iter: the maximum number of iteration
% swarm_size: the population size
% Output:
% Ceq1_fit: the best function value obtained by TLSPSO
% Ceq1: the global best solution obtained by TLSPSO
% valuess: the best value saved with iteration
Ceq1=zeros(1,D);   Ceq1_fit=inf;
Ceq2=zeros(1,D);   Ceq2_fit=inf;
Ceq3=zeros(1,D);   Ceq3_fit=inf;
Ceq4=zeros(1,D);   Ceq4_fit=inf;
NEF=iter*swarm_size;
swarm_size=80;
iter=ceil(NEF/swarm_size);
if size(LB)==1
    LB=repmat(LB,1,D);
end
if size(UB)==1
    UB=repmat(UB,1,D);
end
aa=repmat(LB,swarm_size,1);
bb=repmat(UB,swarm_size,1);
particles_x=aa+rand(swarm_size,D).*(bb-aa);

f_val=feval(fhd,particles_x',varargin{:});
for j=1:swarm_size
    if f_val(j)<Ceq1_fit
        Ceq1_fit=f_val(j);  Ceq1=particles_x(j,:);
        %
    elseif f_val(j)>Ceq1_fit && f_val(j)<Ceq2_fit
        Ceq2_fit=f_val(j);  Ceq2=particles_x(j,:);
    elseif f_val(j)>Ceq1_fit && f_val(j)>Ceq2_fit && f_val(j)<Ceq3_fit
        Ceq3_fit=f_val(j);  Ceq3=particles_x(j,:);
    elseif f_val(j)>Ceq1_fit && f_val(j)>Ceq2_fit && f_val(j)>Ceq3_fit && f_val(j)<Ceq4_fit
        Ceq4_fit=f_val(j);  Ceq4=particles_x(j,:);
        %}
    end
end
p_best=particles_x;
p_best_val=f_val';
g_best_val_t=zeros(1,iter);
valuess=zeros(1,NEF);
lim=zeros(swarm_size,1);
lambda=1;
n1=ceil(swarm_size*0.618);
n2=ceil(swarm_size*1.0);
nef=0;
for i=1:iter
    t=(1-i/iter)^(i/iter);
    fr=0.5*(sin(2*pi*0.25*i)*i/iter+1);
    Ceq_ave=(Ceq1+Ceq2+Ceq3+Ceq4)/4;                              % averaged candidate
    C_pool=[Ceq1; Ceq2; Ceq3; Ceq4; Ceq_ave];
    [~,ind]=sort(p_best_val);
    for j=1:swarm_size
        %
        if (j==ind(n1)) %
            for k=1:D
                rnum=ceil((n1-1)*rand);
                particles_x(j,k)=p_best(ind(rnum),k);
            end
        elseif j==ind(n2) && (i>iter/2)

            n=ceil(5*rand);Ceq=C_pool(n,:);

            if rand<0.015
                particles_x(j,:)=p_best(j,:)+(0.5+0.5*rand)*(Ceq-p_best(j,:));
            else
                particles_x(j,:)=p_best(j,:)+fr*(Ceq-p_best(j,:));
            end
        else
            %}
            r2=rand();                              % r1 and r2 in Eq(15)
            if r2>=t%r2>=0.5
                r=rand(1,D);                                     % r in Eq(11)
                n=ceil(5*rand);
                Ceq=C_pool(n,:);                                  % random selection of one candidate from the pool
                F=2*sign(r-0.5).*(exp(-lambda.*t)-1);%2             % Eq(11)
                r1=rand();
                GCP=0.5*r1*ones(1,D);%0.5                   % Eq(15)
                G0=GCP.*(Ceq-lambda.*p_best(j,:));                      % Eq(14)
                G=G0.*F;                                           % Eq(13)
                particles_x(j,:)=Ceq+(p_best(j,:)-Ceq).*F+(G./lambda).*(1-F);   % Eq(16)
            else
                rnum=ceil(rand*swarm_size);
                rnum1=ceil(rand*swarm_size);while rnum1==rnum,rnum1=ceil(rand*swarm_size);end
                if all(p_best(rnum1,:)==p_best(rnum,:))%
                    for k=1:D
                        if rand<0.5
                            particles_x(j,k)=Ceq1(k);
                        else
                            rnum=ceil((n1-1)*rand);
                            particles_x(j,k)=p_best(ind(rnum),k);

                        end
                    end
                else
                    particles_x(j,:)=p_best(j,:)+(p_best(rnum1,:)-p_best(rnum,:)).*fr;
                end
            end
        end
    end
    particles_x=ControlBoundC(particles_x,aa,bb);
    f_val=feval(fhd,particles_x',varargin{:});
    for j=1:swarm_size
        if f_val(j)<Ceq1_fit
            Ceq1_fit=f_val(j);  Ceq1=particles_x(j,:);            
        elseif f_val(j)>Ceq1_fit && f_val(j)<Ceq2_fit
            Ceq2_fit=f_val(j);  Ceq2=particles_x(j,:);
        elseif f_val(j)>Ceq1_fit && f_val(j)>Ceq2_fit && f_val(j)<Ceq3_fit
            Ceq3_fit=f_val(j);  Ceq3=particles_x(j,:);
        elseif f_val(j)>Ceq1_fit && f_val(j)>Ceq2_fit && f_val(j)>Ceq3_fit && f_val(j)<Ceq4_fit
            Ceq4_fit=f_val(j);  Ceq4=particles_x(j,:);            
        end
        if f_val(j)<p_best_val(j,1)
            p_best(j,:)=particles_x(j,:);
            p_best_val(j,1)=f_val(j);
        else
            lim(j)=lim(j)+1;
        end
    end
    g_best_val_t(i)=Ceq1_fit;
    valuess(nef+1:nef+swarm_size)=Ceq1_fit;
    nef=nef+swarm_size;
    if nef>=NEF
        break
    end
end
valuess=valuess(1:NEF);

function V=ControlBoundC(V,L,U)
%控制变量边界函数
%L为变量低值边界
%U为变量高值边界
%V为待作越界处理的变量
%V,L,U都为相同维数的矩阵
%此函数与V=LimitBound(V,L,U)函数功能相同。
%V=max(V,L);
%V=min(V,U);
c=V<L|V>U;
if sum(c(:))
    V(c)=L(c)+(U(c)-L(c)).*rand(sum(c(:)),1);
end