1.背景

2021年，Braik等人受到变色龙随着环境自适应变化以及寻找食物的动态行为影响，提出了变色龙优化算法(Chameleon Swarm Algorithm，CSA)。

2.算法原理

2.1算法思想

变色龙优化算法模拟了变色龙捕食过程，分为搜索食物-锁定食物-捕获食物三种阶段。

2.2算法过程

群体位置初始化:
$$y=lb+rand\ast (ub-lb)$$
其中，$ub,lb$分别代表变色龙上下位置边界。
搜索食物：
$$y_{t+1}^i=\{\begin{array}{ll}{y}_t^i+p_1(P_t^i-G_t)r_2+p_2(G_t-y_t^i)r_1& r_i⩾Pp\\ y_t^i+\mu ((ub-lb)r_3+lb)\text{sgn}(rand-0.5)& r_i<Pp\end{array}$$
其中，$P_t,G_t$分别变色龙$i$当前最优位置和变色龙群最优位置。

当$r_i\ge Pp$，表示变色龙进行局部探索；当$r_i<pP$进行全局搜索。论文中解释了该公式，假设$P,G,R$是仿射空间（欧式空间去掉了距离）三点，$P,G$两点构成线段中某一点$S$可以由随机数$r_1\in [0,1]$度量：
$$S(r_1)=r_{1}P+(1-r_1)G$$
这里，当$r_1=0$时，此时$S$位于点$G$处；当$r_1=1$时，此时$S$位于点$P$处。同理，点$S,R$中某点$Q$可以由随机数$r_1\in [0,1]$度量：
$$Q(r_1,r_2)=r_{2}S+(1-r_2)R$$
也就是说，通过随机数$r_1,r_2$可以将点$P,R,G$中任意一点表出，这里我们将$Q(r_1,r_2)\leftrightarrow y_{t+1}^i,R\leftrightarrow y_t^i,P\leftrightarrow P_i^t,G\leftrightarrow G_t$就可以将变色龙搜索空间进行表述。
这里，控制参数$\mu$表示变色龙探索能力，这里采用了常用的幂指数控制：
$$\mu =\gamma e^{(-\alpha t/T{)}^{\beta }}$$
随着迭代次数$t$增加，$\mu$不断减少，从而控制算法收敛。
锁定食物:
变色龙通过旋转眼睛锁定食物，论文给出了三维空间中向量平移和旋转推导。

锁定食物时位置更新：
$$y_{t+1}^i=y{r}_t^i\ast {\overline{y}}_t^i$$
这里${\overline{y}}_t^i,y{r}_t^i$分别指变色龙各维度平均值和旋转中心坐标。
$$\{\begin{array}{l}y{r}_t^i=m\times y{c}_t^i\\ y{c}_t^i=y_t^i-y_t^i\\ m=R(\theta ,\overrightarrow{V}{z}_1,\overrightarrow{z_2})\end{array}$$
其中，$y{c}_t^i,m$分别代表定心坐标和旋转矩阵。
捕获食物：
当变色龙靠近食物时，变色龙会进行攻击，速度更新为：
$$v_{t+1}^i=\omega v_t^i+c_1\left(G_t-y_t^i\right)r_1+c_2\left(P_t^i-y_t^i\right)r_2$$
其中，$c_1,c_2$参数控制速度。参数$\omega$更新为：
$$\omega =(1-t/T{)}^{(\rho \sqrt{(}t/T))}$$
位置随着速度变化为：
$$y_{t+1}^i=y_t^i+\left({\left({\nu }_t^i\right)}^2-{\left({\nu }_{t-1}^i\right)}^2\right)/(2a)$$
加速度$a$为：
$$a=2590\times \left(\begin{array}{c}1-e^{-log(t)}\end{array}\right)$$
伪代码

3.代码实现

% 变色龙群算法
function [Best_pos, Best_fitness, Iter_curve, History_pos, History_best]=CSA(pop, dim, ub,lb, fobj, maxIter)
%input
%pop 种群数量
%dim 问题维数
%ub 变量上边界
%lb 变量下边界
%fobj 适应度函数
%maxIter 最大迭代次数
%output
%Best_pos 最优位置
%Best_fitness 最优适应度值
%Iter_curve 每代最优适应度值
%History_pos 每代种群位置
%History_best 每代最优个体位置
%% 初始化种群
for i = 1:dim
    chameleonPositions(:,i) = lb(i)+(ub(i)-lb(i))*rand(pop,1);
end
%% 计算适应度
fit=zeros(pop,1);
for i=1:pop
     fit(i,1)=fobj(chameleonPositions(i,:));
end
Iter_curve=zeros(1,maxIter); %记录每代最优值
fitness=fit; 
[Best_fitness,index]=min(fit);
chameleonBestPosition = chameleonPositions; %初始化变色龙群最优位置
Best_pos = chameleonPositions(index,:); %最优变色龙位置
v=0.1*chameleonBestPosition; %初始速度
v0=0.0*v;
%% 控制参数
rho=1.0;
p1=2.0;  
p2=2.0;  
c1=2.0; 
c2=1.80;  
gamma=2.0; 
alpha = 4.0;  
beta=3.0; 
%% 迭代
for t=1:maxIter
a = 2590*(1-exp(-log(t))); 
omega=(1-(t/maxIter))^(rho*sqrt(t/maxIter)) ; 
p1 = 2* exp(-2*(t/maxIter)^2);  
p2 = 2/(1+exp((-t+maxIter/2)/100)) ;        
mu= gamma*exp(-(alpha*t/maxIter)^beta) ;

% 搜索食物
for i=1:pop
    if rand>=0.1 % pP设置为0.1
        chameleonPositions(i,:)= chameleonPositions(i,:)+ p1*(chameleonBestPosition(i,:)-chameleonPositions(i,:))*rand()+... 
                     + p2*(Best_pos -chameleonPositions(i,:))*rand();
    else 
        chameleonPositions(i,:)= chameleonPositions(i,:)+mu*((ub-lb)*rand+lb)*sign(rand-0.50) ;
    end   
end 

% 捕获食物
     for i=1:pop            
        v(i,:)= omega*v(i,:)+ c1*(chameleonBestPosition(i,:)-chameleonPositions(i,:))*rand +.... 
               + c2*(Best_pos-chameleonPositions(i,:))*rand;        

         chameleonPositions(i,:)=chameleonPositions(i,:)+(v(i,:).^2 - v0(i,:).^2)/(2*a);
     end
    
  v0=v;
  
 %% 边界
 for i=1:pop
     if chameleonPositions(i,:)<lb
        chameleonPositions(i,:)=lb;
     elseif chameleonPositions(i,:)>ub
            chameleonPositions(i,:)=ub;
     end
 end
 
for i=1:pop
    
    ub_=sign(chameleonPositions(i,:)-ub)>0;   
    lb_=sign(chameleonPositions(i,:)-lb)<0;
       
    chameleonPositions(i,:)=(chameleonPositions(i,:).*(~xor(lb_,ub_)))+ub.*ub_+lb.*lb_;  
 
  fit(i,1)=fobj (chameleonPositions(i,:)) ;
      
      if fit(i)<fitness(i)
         chameleonBestPosition(i,:) = chameleonPositions(i,:); 
         fitness(i)=fit(i); 
      end
 end


[fmin,index]=min(fitness); 
if fmin < Best_fitness
    Best_pos = chameleonBestPosition(index,:);
    Best_fitness = fmin;
end

   Iter_curve(t)=Best_fitness; 
   History_pos{t} = chameleonPositions;
   History_best{t} = Best_pos;
end
 
ngPosition=find(fitness== min(fitness)); 
g_best=chameleonBestPosition(ngPosition(1),:);  % Solutin of the problem
Best_fitness =fobj(g_best);

end
%%
function answer = get_orthonormal(m,n)
if ( (nargin==2) && (m>n) && (isnumeric(m)*isnumeric(n)) )
    
elseif ( nargin==1 && isnumeric(m) && length(m)==1 )
    
    n=m;
    
else
   error('Incorrect Inputs. Please read help text in m-file.')
end


count=0;
while (count==0)
    A=rand(m);
    B=A'*A ;


    [P,D] = eig(B) ;

    if ((P'*P - eye(m))>eps) 
        % error, vectors not orthonormal, repeat the random matrix draw again
        count=0;
    else
        % we want the first n of these orthonormal columns
        answer=P(:,1:n) ;
        count=1;
    end


end
end

%%
function [chameleonPositions]=rotation(chameleonPosition, searchAgents, dim)
for i=1:searchAgents      
          if (dim>2) 
              xmax=1;xmin=-1;
              th=round(xmin+rand(1,1)*(xmax-xmin));
              vec=get_orthonormal(dim,2);
              vecA=vec(:,1);
              vecB=vec(:,2);
              theta=(th*rand()*180)*(pi/180) ;
              Rot = RotMatrix(theta,vecA, vecB) ;
             if (theta~=0)
                V=[chameleonPosition(i,:) ]; 
                V_centre=mean(V,1); %Centre, of line
                Vc=V-ones(size(V,1),1)*V_centre; %Centering coordinates

                Vrc=[Rot*Vc']'; %Rotating centred coordinates
%                 Vruc=[Rot*V']'; %Rotating un-centred coordinates
                Vr=Vrc+ones(size(V,1),1)*V_centre; %Shifting back to original location
                 chameleonPosition(i,:)=((Vr)/1); 
 
             end
         else
              xmax=1;xmin=-1;
              th=round(xmin+rand(1,1)*(xmax-xmin));
              theta=th*rand()*180*(pi/180);
              Rot = RotMatrix(theta);
              
               if (theta~=0)
                V=[chameleonPosition(i,:) ];  
                V_centre=mean(V,1); %Centre, of line
                Vc=V-ones(size(V,1),1)*V_centre; %Centering coordinates

                Vrc=[Rot*Vc']'; %Rotating centred coordinates
                Vr=Vrc+ones(size(V,1),1)*V_centre; %Shifting back to original location
                chameleonPosition(i,:)=((Vr)/1);
               end
          end
end  
   chameleonPositions=chameleonPosition;
end
%%
function R = RotMatrix(alpha, u, v)
if numel(alpha) ~= 1
   error('JSimon:RotMatrrix:BadInput1', ...
      'Angle of rotation must be a scalar.');
end

s = sin(alpha);
c = cos(alpha);

% Different algorithms for 2, 3 and N dimensions:
switch nargin
   case 1
      % 2D rotation matrix:
      R = [c, -s;  s, c];
      
   case 2
      if numel(u) ~= 3
         error('JSimon:RotMatrrix:BadAxis2D', ...
            '3D: Rotation axis must have 3 elements.');
      end
      
      % Normalized vector:
      u = u(:);
      u = u ./ sqrt(u.' * u);
            
      % 3D rotation matrix:
      x  = u(1);
      y  = u(2);
      z  = u(3);
      mc = 1 - c;
      R  = [c + x * x * mc,      x * y * mc - z * s,   x * z * mc + y * s; ...
            x * y * mc + z * s,  c + y * y * mc,       y * z * mc - x * s; ...
            x * z * mc - y * s,  y * z * mc + x * s,   c + z * z .* mc];
      
   case 3
      n = numel(u);
      if n ~= numel(v)
         error('JSimon:RotMatrrix:BadAxes3D', ...
            'ND: Axes to define plane of rotation must have the same size.');
      end
      
      % Normalized vectors:
      u = u(:);
      u = u ./ sqrt(u.' * u);
      
      % Care for v being orthogonal to u:
      v = v(:);
      v = v - (u.' * v) * u;
      v = v ./ sqrt(v.' * v);
      
      % Rodrigues' rotation formula:
      R = eye(n) + ...
         (v * u.' - u * v.') * s + ...
         (u * u.' + v * v.') * (c - 1);
      
   otherwise
      error('JSimon:RotMatrrix:BadNInput', ...
            '1 to 3 inputs required.');
end

end

优化问题：
以CEC2005测试函数为例

clear,clc,close all
x = -32:0.1:32;
y = x;
L = length(x);
for i = 1:L
    for j = 1:L
        f(i,j) = fun([x(i) y(j)]);
    end
end
surfc(x, y, f, 'LineStyle', 'none', 'FaceAlpha',0.5);

% 设定变色龙群参数
pop = 50;
dim = 2;
ub = [32, 32];
lb = [-32, -32];
maxIter = 100;
fobj = @(x) fun(x);

% 求解
[Best_pos, Best_fitness, Iter_curve, History_pos, History_best] = CSA(pop, dim, ub, lb, fobj, maxIter);