计算智能——粒子群算法的寻优算法

粒子群算法的寻优算法

1、简介

粒子群优化(PSO, particle swarm optimization)算法是计算智能领域，除了蚁群算法，鱼群算法之外的一种群体智能的优化算法，该算法最早由Kennedy和Eberhart在1995年提出的，该算法源自对鸟类捕食问题的研究。 PSO算法首先在可行解空间中初始化一群粒子，每个粒子都代表极值优化问题的一个潜在最优解，用位置、速度和适应度值三项指标表示该粒子特征。

粒子在解空间中运动，通过跟踪个体极值pbest和群体极值gbest更新个体位置，个体极值pbest是指个体所经历位置中计算得到的适应度值最优位置，群体极值gbest是指种群中的所有粒子搜索到的适应度最优位置。

粒子每更新一次位置，就计算一次适应度值，并且通过比较新粒子的适应度值和个体极值、群体极值的适应度值更新个体极值pbest和群体极值gbest位置。

2、基本原理

在每一次迭代过程中，粒子通过个体极值和群体极值更新自身的速度和位置，更新公式如下:

w为惯性权值，c1,c2为加速系数,pbest为个体最优，gbest为群体最优，pop为当前位置。

基本PSO算法步骤如下:
(1)粒子群初始化;
(2)根据目标函数计算各粒子适应度值，并初始化个体、全局最优值
(3)判断是否满足终止条件，是则搜索停止，输出搜索结果;否则继续下步;
(4)根据速度、位置更新公式更新各粒子的速度和位置;
(5)根据目标函数计算各粒子适应度值;
(6)更新各粒子历史最优值以及全局最优值;(7)跳转至步骤3。
对于终止条件，通常可以设置为适应值误差达到预设要求，或迭代次数超过最大允许迭代次数。
流程图：

群体极值是好找的。但个体极值：第一步时，每个个体的值都是个体极值，第二步后才开始有真正的个体极值的概念。

3、代码实现

1、PSO.m
**设定w函数：w=ws- (ws-we)*(i/maxgen);线性变化惯性权重w（ws,we分别为开始和结束时候的惯性权重）

%% 参数初始化
%粒子群算法中的三个参数
c1 = 1.49445;%加速因子
c2 = 1.49445;
%%c1s=2;
%%c1e=0.5;
%%c2s=0.5;
%%c2e=2;
w = 0.8   %惯性权重
%%ws = 0.9
%%we = 0.4 

maxgen=1000;   % 进化次s数  
sizepop=200;   %种群规模

Vmax=1;       %限制速度围
Vmin=-1;     
popmax=5;    %变量取值范围
popmin=-5;
dim=10;       %适应度函数维数

func=1;       %选择待优化的函数，1为Rastrigin,2为Schaffer,3为Griewank
Drawfunc(func);%画出待优化的函数，只画出二维情况作为可视化输出

%% 产生初始粒子和速度
for i=1:sizepop
    %随机产生一个种群
    pop(i,:)=popmax*rands(1,dim);    %初始种群
    V(i,:)=Vmax*rands(1,dim);             %初始化速度
                                     %计算适应度
    fitness(i)=fun(pop(i,:),func);   %粒子的适应度
end

%% 个体极值和群体极值
[bestfitness bestindex]=min(fitness);
gbest=pop(bestindex,:);   %全局最佳
pbest=pop;                %个体最佳
fitnesspbest=fitness;     %个体最佳适应度值
fitnessgbest=bestfitness; %全局最佳适应度值

%% 迭代寻优
for i=1:maxgen
    
    fprintf('第%d代，',i);
    fprintf('最优适应度%f\n',fitnessgbest);
            w=ws-(ws-we)*(i/maxgen);
            c1=c1s-(c1s-c1e)*(i/maxgen);
            c2=c2s-(c2s-c2e)*(i/maxgen);
    for j=1:sizepop
        
        %速度更新
        V(j,:) = w*V(j,:) + c1*rand*(pbest(j,:) - pop(j,:)) + c2*rand*(gbest - pop(j,:)); %根据个体最优pbest和群体最优gbest计算下一时刻速度
        V(j,find(V(j,:)>Vmax))=Vmax;   %限制速度不能太大
        V(j,find(V(j,:)<Vmin))=Vmin;
        
        %种群更新
        pop(j,:)=pop(j,:)+0.5*V(j,:);       %位置更新
        pop(j,find(pop(j,:)>popmax))=popmax;%坐标不能超出范围
        pop(j,find(pop(j,:)<popmin))=popmin;

        if rand>0.98                         %加入变异种子，用于跳出局部最优值
            pop(j,:)=rands(1,dim);
        end

        %更新第j个粒子的适应度值
        fitness(j)=fun(pop(j,:),func); 
   
    end
    
    for j=1:sizepop
        
        %个体最优更新
        if fitness(j) < fitnesspbest(j)
            pbest(j,:) = pop(j,:);
            fitnesspbest(j) = fitness(j);
        end
        
        %群体最优更新
        if fitness(j) < fitnessgbest
            gbest = pop(j,:);
            fitnessgbest = fitness(j);
        end
    end 
    yy(i)=fitnessgbest;    
        
end
%% 结果分析
figure;
plot(yy)
title('线性递减寻优过程','fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);

2、Rastrigin.m

function y = Rastrigin(x)
% Rastrigin函数
% 输入x,给出相应的y值,在x = ( 0 , 0 ,…, 0 )处有全局极小点0.
[row,col] = size(x);
if  row > 1 
    error( ' 输入的参数错误 ' );
end
y =sum(x.^2-10*cos(2*pi*x)+10);
%y =-y;

3、Schaffer.m

function y=Schaffer(x)

[row,col]=size(x);
if row&g