MATLAB机器学习系列-11：粒子群优化原理及其matlab实现

水深无声 2022-11-21 14:47 365阅读 0赞

## 粒子群优化原理 ##

粒子群算法思想来源于实际生活中鸟捕食的过程。假设在一个n维的空间中，有一群鸟（m只）在捕食，食物位于n维空间的某个点上，对于第i只鸟某一时刻来说，有两个向量描述，一个是鸟的位置向量，第二个是鸟的速度。假设鸟能够判断一个位置的好坏，所谓“好坏”，就是离食物更近了还是更远了。鸟在捕食的过程中会根据自己的经验以及鸟群中的其他鸟的位置决定自己的速度，根据当前的位置和速度，可以得到下一刻的位置，这样每只鸟通过向自己和鸟群学习不断的更新自己的速度位置，最终找到食物，或者离食物足够近的点。更新速度和位置的表达式如下。  
![2020110209224477.png_pic_center][]

![20201102111504923.gif_pic_center][]

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

![20201102093849117.png_pic_center][]

粒子群优化和遗传算法对比  
相同点

*  种群随机初始化
 *  适应度函数值与目标最优解之间的映射

不同点

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

## 案例一：一元函数优化 ##

适应度函数  
**fun.m**

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

**main.m 部分**

**第一步：**

%% I. 清空环境
    clc
    clear all

**第二步：**

x的取值范围1-2

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

![20201102104950277.png_pic_center][]

**第三步：**

%% III. 参数初始化
    c1 = 1.49445;
    c2 = 1.49445;
    
    maxgen = 50;   % 进化次数  
    sizepop = 10;   %种群规模
    
    Vmax = 0.5;%最大速度
    Vmin = -0.5;%最小速度
    popmax = 2;%x边界
    popmin = 1;%x边界

**第四步：**

%% IV. 产生初始粒子和速度
    for i = 1:sizepop %遍历每一个个体
        % 随机产生一个种群
        pop(i,:) = (rands(1) + 1) / 2 + 1;    %初始种群%变化范围是1-2即x的取值
        V(i,:) = 0.5 * rands(1);  %初始化速度
        % 计算适应度
        fitness(i) = fun(pop(i,:));   
    end

种子pop  
![20201102102550191.png_pic_center][]

初始化速度  
![20201102102818529.png_pic_center][]

适应度函数值  
![20201102102702461.png_pic_center][]

**第五步：**

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

在本例子中bestindex=1,因为是随机化，所以每次运行不一样

![20201102103524339.png_pic_center][] ![20201102103450396.png_pic_center][] ![2020110210330677.png_pic_center][] ![20201102103414371.png_pic_center][]

**第六步：**

%% 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;%X的最大边界限制
            pop(j,find(pop(j,:)<popmin)) = popmin;%x的最小边界限制
            
            % 适应度值更新
            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;   %一轮的群体最优适应度函数值赋值给yy(i)       
    end

yy有50个值

![20201102104737255.png_pic_center][] 我们从图中可以发现fitnesszbest，zbest 是最优值（对比前面出现的函数图） ![20201102104909333.png_pic_center][]

第七步：

%% VII. 输出结果并绘图
    [fitnesszbest zbest]
    plot(zbest, fitnesszbest,'r*')
    
    figure
    plot(yy)
    title('最优个体适应度','fontsize',12);
    xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);

从图中发现到20轮时就基本找到最大值啦。  
![20201102105253274.png_pic_center][]

![20201102105337425.png_pic_center][]

## 案例二：二元函数优化 ##

二元函数代码和前面一元函数类似  
fun.m

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;

main.m

%% 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);

![20201102110820130.png_pic_center][] ![20201102110847530.png_pic_center][]

MATLAB 2014b之后好像有粒子群的工具箱。自己去研究吧。  
MATLAB非我主语言。  
[MATLAB-pso工具箱][MATLAB-pso]  
![在这里插入图片描述][20201102110911696.jpg_pic_center]

作者：电气工程的计算机萌新-余登武  
![在这里插入图片描述][20201102110953211.gif_pic_center]

[2020110209224477.png_pic_center]: /images/20221120/da76364b6d744725a366218b007518e2.png
[20201102111504923.gif_pic_center]: /images/20221120/de3c52cc13fc48f096917378601cecb1.png
[20201102093849117.png_pic_center]: /images/20221120/3e36f2f504624eb7b3bfbbedac4bc17d.png
[20201102104950277.png_pic_center]: /images/20221120/ec528ae4b817412c970d3c838587aee5.png
[20201102102550191.png_pic_center]: /images/20221120/8f9b573bd55545abb88601fbff7ff7d8.png
[20201102102818529.png_pic_center]: /images/20221120/1164d9d597f8460da3474185ab371a08.png
[20201102102702461.png_pic_center]: /images/20221120/1f2b4ef9a8874c8194b321194f7ad2c9.png
[20201102103524339.png_pic_center]: https://img-blog.csdnimg.cn/20201102103524339.png#pic_center
[20201102103450396.png_pic_center]: https://img-blog.csdnimg.cn/20201102103450396.png#pic_center
[2020110210330677.png_pic_center]: https://img-blog.csdnimg.cn/2020110210330677.png#pic_center
[20201102103414371.png_pic_center]: /images/20221120/2cfc898b4fd443c098a70d8e4f190ff5.png
[20201102104737255.png_pic_center]: https://img-blog.csdnimg.cn/20201102104737255.png#pic_center
[20201102104909333.png_pic_center]: /images/20221120/3afae6d5646f46d4afcbce748952701d.png
[20201102105253274.png_pic_center]: /images/20221120/14012210462f499fbb87a3ca9ae7f99e.png
[20201102105337425.png_pic_center]: /images/20221120/bb34421b524a4ffb8c9351d1fefa599d.png
[20201102110820130.png_pic_center]: https://img-blog.csdnimg.cn/20201102110820130.png#pic_center
[20201102110847530.png_pic_center]: /images/20221120/2784a2da7e4946aea639f82ee9bad2a6.png
[MATLAB-pso]: https://www.mathworks.com/help/gads/particleswarm.html
[20201102110911696.jpg_pic_center]: /images/20221120/0ed6f1a3ba8b43ffabeb7565bbcda8d6.png
[20201102110953211.gif_pic_center]: /images/20221120/2246cadd46cd4d9ca2cabf42b327a230.png