一、麻雀演算法簡介

麻雀搜尋演算法(Sparrow Search Algorithm, SSA)是於2020年提出的。SSA 主要是受麻雀的覓食行為和反捕食行為的啟發而提出的。該演算法比較新穎，具有尋優能力強，收斂速度快的優點。
1 演算法原理
建立麻雀搜尋演算法的數學模型，主要規則如下所述：
（1）發現者通常擁有較高的能源儲備並且在整個種群中負責搜尋到具有豐富食物的區域，為所有的加入者提供覓食的區域和方向。在模型建立中能量儲備的高低取決於麻雀個體所對應的適應度值(Fitness Value)的好壞。
（2）一旦麻雀發現了捕食者，個體開始發出鳴叫作為報警訊號。當報警值大於安全值時，發現者會將加入者帶到其它安全區域進行覓食。
（3）發現者和加入者的身份是動態變化的。只要能夠尋找到更好的食物來源，每隻麻雀都可以成為發現者，但是發現者和加入者所佔整個種群數量的比重是不變的。也就是說，有一隻麻雀變成發現者必然有另一隻麻雀變成加入者。
（4）加入者的能量越低，它們在整個種群中所處的覓食位置就越差。一些飢腸轆轆的加入者更有可能飛往其它地方覓食，以獲得更多的能量。
（5）在覓食過程中，加入者總是能夠搜尋到提供最好食物的發現者，然後從最好的食物中獲取食物或者在該發現者周圍覓食。與此同時，一些加入者為了增加自己的捕食率可能會不斷地監控發現者進而去爭奪食物資源。
（6）當意識到危險時，群體邊緣的麻雀會迅速向安全區域移動，以獲得更好的位置，位於種群中間的麻雀則會隨機走動，以靠近其它麻雀。
在模擬實驗中，我們需要使用虛擬麻雀進行食物的尋找，由n只麻雀組成的種群可表示為如下形式：

其中，d表示待優化問題變數的維數，n則是麻雀的數量。那麼，所有麻雀的適應度值可以表示為如下形式：

其中，f表示適應度值。
在SSA中， 具有較好適應度值的發現者在搜尋過程中會優先獲取食物。此外， 因為發現者負責為整個麻雀種群尋找食物併為所有加入者提供覓食的方向。因此，發現者可以獲得比加入者更大的覓食搜尋範圍。根據規則(1)和規則(2)，在每次迭代的過程中，發現者的位置更新描述如下：

其中， t代表當前迭代數， j=1， 2， 3， ...， d.item maz是一個常數，表示最大的迭代次數。Xy表示第i個麻雀在第j維中的位置資訊。xE(0，1]是一個隨機數。R2(R2E[0，1])和ST(STe[0.5， 1] ) 分別表示預警值和安全值。Q是服從正態分佈的隨機數.L表示一個1×d的矩陣， 其中該矩陣內每個元素全部為
1.
當R2<ST時，這意味著此時的覓食環境周圍沒有捕食者，發現者可以執行廣泛的搜尋操作。如果R2≥ST，這表示種群中的一些麻雀已經發現了捕食者，井向種群中其它麻雀發出了警報，此時所有麻雀都需要迅速飛到其它安全的地方進行覓食。對於加入者，它們需要執行規則(3)和規則(4)。如前面所描述，在覓食過程中，一些加入者會時刻監視著發現者。一旦它們察覺到發現者已經找到了更好的食物，它們會立即離開現在的位置去爭奪食物。如果它們贏了，它們可以立即獲得該發現者的食物，否則需要繼續執行規則(4)。加入者的位置更新描述如下：

其中， X， 是目前發現者所佔據的最優位置， X worst則表示當前全域性最差的位置。A表示一個1×d的矩陣， 其中每個元素隨機賦值為1或-1，並且A+=A(AA)-.當i>n/2時，這表明，適應度值較低的第i個加入者沒有獲得食物，處於十分飢餓的狀態，此時需要飛往其它地方覓食，以獲得更多的能量。在模擬實驗中，我們假設這些意識到危險的麻雀佔總數量的10%到20%。這些麻雀的初始位置是在種群中隨機產生的。根據規則(5)，其數學表示式可以表示為如下形式：

其中， 其中X best是當前的全域性最優位置。β作為步長控制引數， 是服從均值為0， 方差為1的正態分佈的隨機數。KE[-1， 1] 是一個隨機數，則是當前麻雀個體的適應度值。f，和fw分別是當前全域性最佳和最差的適應度值。e的常數，以避免分母出現零。
為簡單起見， 當f：>f， 表示此時的麻雀正處於種群的邊緣， 極其容易受到捕食者的攻擊。X best表示這個位置的麻雀是種群中最好的位置也是十分安全的。f；=f，時，這表明處於種群中間的麻雀意識到了危險，需要靠近其它的麻雀以此儘量減少它們被捕食的風險。K表示麻雀移動的方向同時也是步長控制引數。

2 演算法流程

Step1： 初始化種群，迭代次數，初始化捕食者和加入者比列。
Step2：計算適應度值，並排序。
Step3：利用式（3）更新捕食者位置。
Step4：利用式（4）更新加入者位置。
Step5：利用式（5）更新警戒者位置。
Step6：計算適應度值並更新麻雀位置。
Step7：是否滿足停止條件，滿足則退出，輸出結果，否則，重複執行Step2-6；

二、原始碼

% -----------------------------------------------------------------------------------------------------------
% Sparrow Search algorithm (SSA) (demo)



clear all 
clc

SearchAgents_no=100; % 代數

Function_name='F1'; % Name of the test function that can be from F1 to F23 (Table 1,2,3 in the paper)

Max_iteration=1000; % Maximum numbef of iterations

% Load details of the selected benchmark function
[lb,ub,dim,fobj]=Get_Functions_details(Function_name);

[fMin,bestX,SSA_curve]=SSA(SearchAgents_no,Max_iteration,lb,ub,dim,fobj);  


%Draw objective space
 semilogy(SSA_curve,'Color','g')
  
axis ([0 1000 0 1 ])
title('Objective space')
xlabel('Iteration');
ylabel('Best score obtained so far');
%axis tight
grid on
box on
legend('SSA')
display(['The best solution obtained by SSA is : ', num2str(bestX)]);
  display(['The best optimal value of the objective funciton found by SSA is : ', num2str(fMin)]);

%pop是種群，M是迭代次數，fobj是用來計算適應度的函式
%pNum是生產者
%r1==a
% r2:預警值
% fit：適應度

function [fMin , bestX,Convergence_curve ] = SSA(pop, M,c,d,dim,fobj  )
        
   P_percent = 0.2;    % The population size of producers accounts for "P_percent" percent of the total population size       


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pNum = round( pop *  P_percent );    % The population size of the producers   


lb= c.*ones( 1,dim );    % Lower limit/bounds/     a vector
ub= d.*ones( 1,dim );    % Upper limit/bounds/     a vector
%Initialization
for i = 1 : pop
    
    x( i, : ) = lb + (ub - lb) .* rand( 1, dim );  
    fit( i ) = fobj( x( i, : ) ) ;                       
end
pFit = fit;                      
pX = x;                            % The individual's best position corresponding to the pFit
[ fMin, bestI ] = min( fit );      % fMin denotes the global optimum fitness value
bestX = x( bestI, : );             % bestX denotes the global optimum position corresponding to fMin
 

 % Start updating the solutions.
%
for t = 1 : M    
  
      
  [ ans, sortIndex ] = sort( pFit );% Sort.
     
  [fmax,B]=max( pFit );
   worse= x(B,:);  
         
   r2=rand(1);
   
    %%%%%%%%%%%%%5%%%%%%這一部位為發現者（探索者）的位置更新%%%%%%%%%%%%%%%%%%%%%%%%%  
if(r2<0.8)%預警值較小，說明沒有捕食者出現
 
    for i = 1 : pNum  %r2小於0.8的發現者的改變（1-20）                                                 % Equation (3)
         r1=rand(1);
        x( sortIndex( i ), : ) = pX( sortIndex( i ), : )*exp(-(i)/(r1*M));%對自變數做一個隨機變換
        x( sortIndex( i ), : ) = Bounds( x( sortIndex( i ), : ), lb, ub );%對超過邊界的變數進行去除
        fit( sortIndex( i ) ) = fobj( x( sortIndex( i ), : ) );   %就算新的適應度值
    end
  else   %預警值較大，說明有捕食者出現威脅到了種群的安全，需要去其它地方覓食
  for i = 1 : pNum   %r2大於0.8的發現者的改變
          
  x( sortIndex( i ), : ) = pX( sortIndex( i ), : )+randn(1)*ones(1,dim);
  x( sortIndex( i ), : ) = Bounds( x( sortIndex( i ), : ), lb, ub );
  fit( sortIndex( i ) ) = fobj( x( sortIndex( i ), : ) );
       
  end
      
end
 
 
        

      %%%%%%%%%%%%%5%%%%%%這一部位為加入者（追隨者）的位置更新%%%%%%%%%%%%%%%%%%%%%%%%%
   for i = ( pNum + 1 ) : pop     %剩下20-100的個體的變換                % Equation (4)
     
         A=floor(rand(1,dim)*2)*2-1;
         
          if( i>(pop/2))%這個代表這部分麻雀處於十分飢餓的狀態（因為它們的能量很低，也是是適應度值很差），需要到其它地方覓食
           x( sortIndex(i ), : )=randn(1)*exp((worse-pX( sortIndex( i ), : ))/(i)^2);
          else%這一部分追隨者是圍繞最好的發現者周圍進行覓食，其間也有可能發生食物的爭奪，使其自己變成生產者
        x( sortIndex( i ), : )=bestXX+(abs(( pX( sortIndex( i ), : )-bestXX)))*(A'*(A*A')^(-1))*ones(1,dim);  

         end  
        x( sortIndex( i ), : ) = Bounds( x( sortIndex( i ), : ), lb, ub );%判斷邊界是否超出
        fit( sortIndex( i ) ) = fobj( x( sortIndex( i ), : ) );%計算適應度值
        
   end
  %%%%%%%%%%%%%5%%%%%%這一部位為意識到危險（注意這裡只是意識到了危險，不代表出現了真正的捕食者）的麻雀的位置更新%%%%%%%%%%%%%%%%%%%%%%%%%
  c=randperm(numel(sortIndex));%%%%%%%%%這個的作用是在種群中隨機產生其位置（也就是這部分的麻雀位置一開始是隨機的，意識到危險了要進行位置移動，
                                                                         %處於種群外圍的麻雀向安全區域靠攏，處在種群中心的麻雀則隨機行走以靠近別的麻雀）


        x( sortIndex( b(j) ), : )=bestX+(randn(1,dim)).*(abs(( pX( sortIndex( b(j) ), : ) -bestX)));

    else                       %處於種群中心的麻雀的位置改變

        x( sortIndex( b(j) ), : ) =pX( sortIndex( b(j) ), : )+(2*rand(1)-1)*(abs(pX( sortIndex( b(j) ), : )-worse))/ ( pFit( sortIndex( b(j) ) )-fmax+1e-50);

          end
        x( sortIndex(b(j) ), : ) = Bounds( x( sortIndex(b(j) ), : ), lb, ub );
       
       fit( sortIndex( b(j) ) ) = fobj( x( sortIndex( b(j) ), : ) );
 end
    for i = 1 : pop 
        if ( fit( i ) < pFit( i ) )
            pFit( i ) = fit( i );
            pX( i, : ) = x( i, : );
        end
        

            
        end
    end
  
    Convergence_curve(t)=fMin;
  
end


% Application of simple limits/bounds
function s = Bounds( s, Lb, Ub)
  % Apply the lower bound vector

  
  % Apply the upper bound vector 
  J = temp > Ub;
  temp(J) = Ub(J);
  % Update this new move 
  s = temp;

%---------------------------------------------------------------------------------------------------------------------------



 

三、執行結果

四、matlab版本及參考文獻

1 matlab版本
2014a

2 參考文獻
[1] 包子陽,餘繼周,楊杉.智慧優化演算法及其MATLAB例項（第2版）[M].電子工業出版社，2016.
[2]張巖,吳水根.MATLAB優化演算法原始碼[M].清華大學出版社，2017.
[3]周品.MATLAB 神經網路設計與應用[M].清華大學出版社，2013.
[4]陳明.MATLAB神經網路原理與例項精解[M].清華大學出版社，2013.
[5]方清城.MATLAB R2016a神經網路設計與應用28個案例分析[M].清華大學出版社，2018.