程序代码：

popsize=20;                 %群体规模
MAXITER=1000;               %最大迭代次数
dimension=2;               %问题维数
irange_l=-100;              %位置初始化下界
irange_r=100;               %位置初始化上界
xmax=100;                   %搜索范围上界
xmin=0;                     %搜索范围下界
M=(xmax-xmin)/2;            %搜索范围中值
sum1=0;
runno=50;                   %算法运行次数
data1=zeros(runno,MAXITER); %记录每一轮中每一步迭代的最好适应值
bestVariable=zeros(runno,dimension);   %记录每次运行时取得最优适应值的变量值
T=cputime;              %记录CPU时间
 
for run=1:runno             %循环    
    x=(irange_r-irange_l)*rand(popsize,dimension,1)+irange_l;     %初始化粒子当前位置
    pbest=x;                %将粒子个体最好位置初始化为当前最好位置
    gbest=zeros(1,dimension);   %初始化全局最好位置变量
    f_pbest=zeros(1,popsize);   %粒子个体最好适应值
    f_x=zeros(1,popsize);       %粒子个体位置适应值
     
    for i=1:popsize             %计算当前位置和个体最好位置的适应值
        f_x(i)=F(x(i,:));
        f_pbest(i)=f_x(i);
    end    
    g=min(find(f_pbest==min(f_pbest(1:popsize)),1));  %找到全局最好位置的粒子下标
    gbest=pbest(g,:);                               %找到全局最好位置
    f_gbest=f_pbest(g);                             %记录全局最好位置的适应值
    MINIMUM=f_gbest;                             %将全局最好位置的适应值记录为MINIMUM
     %初始化最佳变量值
    for dim=1:dimension
         bestVariable(run,dim)=x(g,dim);
    end
   
    for t=1:MAXITER                                 %每轮计算中，算法的迭代过程
        alpha=(1.0-0.5)*(MAXITER-t)/MAXITER+0.5;    %收缩-扩展系数
        mbest=sum(pbest)/popsize;                   %计算平均最好位置
        for i=1:popsize
            fi1=rand(1,dimension);
            fi2=rand(1,dimension);
            p=(2*fi1.*pbest(i,:)+2.1*fi2.*gbest)./(2*fi1+2.1*fi2);    %计算随机点
            u=rand(1,dimension);
            b=alpha*abs(mbest-x(i,:));
            v=-log(u);
            x(i,:)=p+((-1).^ceil(0.5+rand(1,dimension))).*b.*v;     %粒子位置的更新
             
            %%%%以下是将粒子位置限制在收缩范围之内
            z=x(i,:)-(xmin+xmax)/2;
            z=sign(z).*min(abs(z),M);
            x(i,:)=z+(xmax+xmin)/2;  
             
            f_x(i)=F(x(i,:));                        %计算粒子当前位置适应值
            if f_x(i)<f_pbest(i)                     
                pbest(i,:)=x(i,:);                   %更新粒子个体最好位置
                f_pbest(i)=f_x(i);                   %更新粒子个体最好适应值
            end
             
            if f_pbest(i)<f_gbest                    
                gbest=pbest(i,:);                    %更新粒子群全局最好位置
                f_gbest=f_pbest(i);                  %更新粒子群全局最优适应值
                %记录每次运行的最佳变量值
                for dim=1:dimension
                    bestVariable(run,dim)=x(i,dim);
                end                
            end
            MINIUM=f_gbest;                          %记录全局最好位置适应值       
        end
        data1(run,t)=MINIUM;                         %记录本次迭代找到的最好适应值
    end    
    sum1=sum1+MINIUM;   
end
st=cputime-T;                                      %计算所有计算所需时间
aveRunTime=st/runno;                               %每轮使用的平均CPU时间
aveValue=sum1/runno;                               %计算runno论得到的最好适应值的平均值
% x=1:MAXITER;
% plot(x,data1(1,:))