粒子群优化算法PSO及matlab实现

算法学习自：MATLAB与机器学习教学视频

1、粒子群优化算法概述

粒子群优化（PSO, particle swarm optimization）算法是计算智能领域，除了蚁群算法，鱼群算法之外的一种群体智能的优化算法，该算法最早由Kennedy和Eberhart在1995年提出的，该算法源自对鸟类捕食问题的研究。
　　• PSO算法首先在可行解空间中初始化一群粒子，每个粒子都代表极值优化问题的一个潜在最优解，用位置、速度和适应度值三项指标表示该粒子特征。
　　• 粒子在解空间中运动，通过跟踪个体极值Pbest和群体极值Gbest更新个体位置，个体极值Pbest是指个体所经历位置中计 算得到的适应度值最优位置，群体极值Gbest是指种群中的所有粒子搜索到的适应度最优位置。
　　• 粒子每更新一次位置，就计算一次适应度值，并且通过比较新粒子的适应度值和个体极值、群体极值的适应度值更新个体 极值Pbest和群体极值Gbest位置。

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

　　　　　　　

　　　　　　

 2、粒子群优化算法与遗传算法对比

• 相同点：
　　种群随机初始化
　　适应度函数值与目标最优解之间的映射
• 不同点：
　　PSO算法没有选择、交叉、变异等操作算子
　　PSO有记忆的功能
　　信息共享机制不同，遗传算法是互相共享信息，整个种群的移动是比较均匀地向最优区域移动，而在PSO中，只 有gBest或lBest给出信息给其他粒子，属于单向的信息流动，整个搜索更新过程是跟随当前最优解的过程。因此， 在一般情况下，PSO的收敛速度更快。

3、案例分析

①一元函数优化

　　在该例中，想要优化的一元函数为，当x范围在[1,2]，寻找它的极大值。

matlab代码如下：

%% I. 清空环境
clc
clear all

%% II. 绘制目标函数曲线图
x = 1:0.01:2;
y = sin(10*pi*x) ./ x;
figure
plot(x, y)
hold on

%% III. 参数初始化
c1 = 1.49445;
c2 = 1.49445;

maxgen = 50;   % 进化次数  
sizepop = 10;   %种群规模

Vmax = 0.5;
Vmin = -0.5;
popmax = 2;
popmin = 1;

%% IV. 产生初始粒子和速度
for i = 1:sizepop
    % 随机产生一个种群
    pop(i,:) = (rands(1) + 1) / 2 + 1;    %初始种群
    V(i,:) = 0.5 * rands(1);  %初始化速度
    % 计算适应度
    fitness(i) = fun(pop(i,:));   
end

%% V. 个体极值和群体极值
[bestfitness, bestindex] = max(fitness);
zbest = pop(bestindex,:);   %全局最佳
gbest = pop;    %个体最佳
fitnessgbest = fitness;   %个体最佳适应度值
fitnesszbest = bestfitness;   %全局最佳适应度值

%% VI. 迭代寻优
for i = 1:maxgen
    
    for j = 1:sizepop
        % 速度更新
        V(j,:) = V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));
        V(j,V(j,:)>Vmax) = Vmax;
        V(j,V(j,:)<Vmin) = Vmin;
        
        % 种群更新
        pop(j,:) = pop(j,:) + V(j,:);
        pop(j,pop(j,:)>popmax) = popmax;
        pop(j,pop(j,:)<popmin) = popmin;
        
        % 适应度值更新
        fitness(j) = fun(pop(j,:)); 
    end
    
    for j = 1:sizepop    
        % 个体最优更新
        if fitness(j) > fitnessgbest(j)
            gbest(j,:) = pop(j,:);
            fitnessgbest(j) = fitness(j);
        end
        
        % 群体最优更新
        if fitness(j) > fitnesszbest
            zbest = pop(j,:);
            fitnesszbest = fitness(j);
        end
    end 
    yy(i) = fitnesszbest;          
end

%% VII. 输出结果并绘图
[fitnesszbest zbest];
plot(zbest, fitnesszbest,'r*')

figure
plot(yy)
title('最优个体适应度','fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);

function y = fun(x)
% 函数用于计算粒子适应度值
%x           input           输入粒子 
%y           output          粒子适应度值 
y = sin(10 * pi * x) / x;

结果图示：

2、二元函数优化

 　　在该例中，想要优化的二元函数为，当x和y范围都在[-5,5]，寻找它的极大值。

matlab代码如下：

function y = fun(x)
%函数用于计算粒子适应度值
%x           input           输入粒子 
%y           output          粒子适应度值 
y = x(1).^2 + x(2).^2 - 10*cos(2*pi*x(1)) - 10*cos(2*pi*x(2)) + 20;

%% I. 清空环境
clc
clear

%% II. 绘制目标函数曲线
figure
[x,y] = meshgrid(-5:0.1:5,-5:0.1:5);
z = x.^2 + y.^2 - 10*cos(2*pi*x) - 10*cos(2*pi*y) + 20;
mesh(x,y,z)
hold on

%% III. 参数初始化
c1 = 1.49445;
c2 = 1.49445;

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

Vmax = 1;
Vmin = -1;
popmax = 5;
popmin = -5;

%% IV. 产生初始粒子和速度
for i = 1:sizepop
    % 随机产生一个种群
    pop(i,:) = 5*rands(1,2);    %初始种群
    V(i,:) = rands(1,2);  %初始化速度
    % 计算适应度
    fitness(i) = fun(pop(i,:));   %染色体的适应度
end

%% V. 个体极值和群体极值
[bestfitness bestindex] = max(fitness);
zbest = pop(bestindex,:);   %全局最佳
gbest = pop;    %个体最佳
fitnessgbest = fitness;   %个体最佳适应度值
fitnesszbest = bestfitness;   %全局最佳适应度值

%% VI. 迭代寻优
for i = 1:maxgen
    
    for j = 1:sizepop
        % 速度更新
        V(j,:) = V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));
        V(j,find(V(j,:)>Vmax)) = Vmax;
        V(j,find(V(j,:)<Vmin)) = Vmin;
        
        % 种群更新
        pop(j,:) = pop(j,:) + V(j,:);
        pop(j,find(pop(j,:)>popmax)) = popmax;
        pop(j,find(pop(j,:)<popmin)) = popmin;
        
        % 适应度值更新
        fitness(j) = fun(pop(j,:)); 
    end
    
    for j = 1:sizepop  
        % 个体最优更新
        if fitness(j) > fitnessgbest(j)
            gbest(j,:) = pop(j,:);
            fitnessgbest(j) = fitness(j);
        end
        
        % 群体最优更新
        if fitness(j) > fitnesszbest
            zbest = pop(j,:);
            fitnesszbest = fitness(j);
        end
    end 
    yy(i) = fitnesszbest;            
end
%% VII.输出结果
[fitnesszbest, zbest]
plot3(zbest(1), zbest(2), fitnesszbest,'bo','linewidth',1.5)

figure
plot(yy)
title('最优个体适应度','fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);

结果图示：

    

4、速度更新权重W的选择

示例代码如下：

ws = 0.9;
we = 0.4;
maxgen = 300;
hold on;

for k = 1:maxgen
    w(k) = ws - (ws-we)*(k/maxgen);
end
plot(w,'linewidth',1.5);

for k = 1:maxgen
    w(k) = ws - (ws-we)*(k/maxgen)^2;
end
plot(w,'r-.','linewidth',1.5);

for k = 1:maxgen
    w(k) = ws - (ws-we)*(2*k/maxgen-(k/maxgen)^2);
end
plot(w,'g:','linewidth',1.5);

for k = 1:maxgen
    w(k) = we * (ws/we)^(1/(1+10*k/maxgen));
end
plot(w,'y--','linewidth',1.5);

legend('Rule-1','Rule-2','Rule-3','Rule-4')
xlabel('迭代次数')
ylabel('速度更新权重W')

图示：