粒子群算法简介

最新推荐文章于 2024-07-16 14:14:36 发布

Robert_YHQ

最新推荐文章于 2024-07-16 14:14:36 发布

阅读量1.6k

点赞数 1

文章标签：算法 matlab

本文链接：https://blog.csdn.net/Robert_YHQ/article/details/121216934

版权

粒子群优化(Particle Swarm Optimization, PSO)算法是Kennedy和Eberhart受人工生命研究结果的启发、通过模拟鸟群觅食过程中的迁徙和群聚行为而提出的一种基于群体智能的全局随机搜索算法，自然界中各种生物体均具有一定的群体行为，而人工生命的主要研究领域之一是探索自然界生物的群体行为，从而在计算机上构建其群体模型。自然界中的鸟群和鱼群的群体行为一直是科学家的研究兴趣，生物学家Craig Reynolds在1987年提出了一个非常有影响的鸟群聚集模型，在他的仿真中，每一个个体遵循：

(1) 避免与邻域个体相冲撞；

(2) 匹配邻域个体的速度；

(3) 飞向鸟群中心，且整个群体飞向目标。

仿真中仅利用上面三条简单的规则，就可以非常接近的模拟出鸟群飞行的现象。1995年，美国社会心理学家James Kennedy和电气工程师Russell Eberhart共同提出了粒子群算法，其基本思想是受对鸟类群体行为进行建模与仿真的研究结果的启发。他们的模型和仿真算法主要对Frank Heppner的模型进行了修正，以使粒子飞向解空间并在最好解处降落。

相比于其他算法，PSO有三个比较显著的优点：1.简单易行。2.收敛速度快。3.设置参数少。

PSO应用的问题：多目标优化、分类、模式识别、决策。

每个寻优的问题都被看成是一个“粒子”，所有粒子都在一个空间中进行搜索最优解。所有的粒子都由一个函数来判断当前的位置的好坏，并且每一个粒子都应赋予储存的功能，以便能得到最优解，每一个粒子都需要有自己的速度，速度是由粒子本身和群体来进行动态调整的。故我们在D维空间中建立D维向量，则每一个粒子可以用一个向量来表示其位置，记作 $x_{i}=\left ( x_{i1},x_{i2},x_{i3},...x_{iD} \right )$

同样地，还会有一个向量来表示其速度，记作 $v_{i}=\left ( v_{i1},v_{i2},v_{i3},...v_{iD} \right )$ ,每个个体会经过一个最好的位置，记作 $pbest_{i}$ ，同时，粒子群体会经过一个最好的位置，记作 $gbest_{}$ 。每个粒子的位置和速度都会被限制在一定的范围之内。

粒子i的第d维速度更新公式：

$v_{id}^{k}=wv_{id}^{k-1}+c_{1}r_{1}\left ( pbest_{id}-x_{id}^{k-1} \right )+c_{2}r_{2}\left ( gbest_{d} -x_{id}^{k-1}\right )$

粒子i的第d维位置更新公式：

$x_{id}^{k}=x_{id}^{k-1}+v_{id}^{k-1}$

其中，各公式中的k表示当前状态，k-1表示上一状态，类似于数字逻辑中触发器的次态方程的表达形式。 $c_{1},c_{2}$ 表示加速度常量，是调节学习的最大步长。r1，r2表示两个0到1之间的随机数，以增加搜索的随机性，避免因个人因素使结果更偏向于两者之一的结果。w则表示惯性权重，非负数，调节对空间的搜索范围。

速度更新和位置更新受到的影响因素可以很清晰的在公式中表示出来，这里不多说了emmmm。

算法流程：（由于作者能力有限，流程图就算了吧）

1.随机初始化例子群体，包括位置和速度。

2.根据fitness 函数，评价每个粒子的适应度（即当前位置的优劣程度）。

3.根据每个粒子当前的位置来对个体、群体最佳位置进行修正。

4.根据上述公式更新每个粒子的速度和位置。

5.判断是否终止，没终止则返回第二步（终止条件：计数器达到最大迭代次数或全局最优解满足最小界限）

各种构成要素要选择合适的数值进行模拟，过大过小的话emmmmm，懂得都懂。

这里给惯性权重单独拎出来，惩罚它：

惯性权重w描述粒子上一状态速度对挡墙速度的影响。当w值较大时，全局寻优能力强，局部寻优能力较弱；反之则局部寻优能力强。当问题空间较大时，为了在搜索速度和搜索精度之间达到平衡，通常做法是使算法在前期有较高的全局搜索能力，在后期具有较高的局部搜索能力，故w不适合是一个固定的常数。公式如下：

$w=w_{max}-\left ( w_{max}-w_{min} \right )*\frac{run}{run_{max}}$

顾名思义max和min代表着最大和最小，run代表着当前迭代次数。

其他参数则需自我设置。

代码：（网上找的+自己改了一部分）

%------初始格式化--------------------------------------------------
 
clear all;
clc;
format long;
%------给定初始化条件----------------------------------------------
 
c1=2;             %学习因子1
c2=2;             %学习因子2
w=0.7298;              %惯性权重
MaxDT=200;            %最大迭代次数
% D=2;                  %搜索空间维数（未知数个数）
 
N=20;                  %初始化群体个体数目
%eps=10^(-6);           %设置精度(在已知最小值时候用)
 
Vmax=1;
Vmin=-1;
popmax=5;
popmin=-5;
%------初始化种群的个体(可以在这里限定位置和速度的范围)------------
 
for i=1:N
        pop(i,:)=popmin+(popmax-popmin)*rand(1,2);  %随机初始化位置
 
        V(i,:)=rand(1,2); %随机初始化速度
        fitness(i)=ackley(pop(i,:));
 
    end
%------先计算各个粒子的适应度，并初始化Pi和Pg----------------------
 
[fitnessgbest bestindex]=min(fitness);
gbest=pop(bestindex,:);
pbest=pop;
fitnesspbest=fitness;
 
for i=1:MaxDT
    for j=1:N
        V(j,:)=w*V(j,:)+c1*rand*(pbest(j,:)-pop(j,:))+c2*rand*(gbest-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;
        
%         if rand>0.8
 
%             k=ceil(2*rand);
 
%             pop(j,k)=rand;
 
%         end        
 
        fitness(j)=ackley(pop(j,:));
  
    
       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
%------最后给出计算结果
 
plot(yy)
title(['适应度曲线 ' '终止次数=' num2str(MaxDT)]);
xlabel('进化代数');
ylabel('适应度')
%------算法结束---DreamSun GL & HF-----------------------------------

优化函数：

% ackley.m
 
% Ackley's function, from http://www.cs.vu.nl/~gusz/ecbook/slides/16
 
% and further shown at: 
 
% http://clerc.maurice.free.fr/pso/Semi-continuous_challenge/Semi-continuous_challenge.htm
 
%
 
% commonly used to test optimization/global minimization problems
 
%
 
% f(x)= [ 20 + e ...
 
%        -20*exp(-0.2*sqrt((1/n)*sum(x.^2,2))) ...
 
%        -exp((1/n)*sum(cos(2*pi*x),2))];
 
%
 
% dimension n = # of columns of input, x1, x2, ..., xn
 
% each row is processed independently,
 
% you can feed in matrices of timeXdim no prob
 
%
 
% example: cost = ackley([1,2,3;4,5,6])
 
 
function [out]=ackley(in)
 
% dimension is # of columns of input, x1, x2, ..., xn
 
 n=length(in(1,:));
 
 x=in;
 e=exp(1);
 
 out = (20 + e ...
       -20*exp(-0.2*sqrt((1/n).*sum(x.^2,2))) ...
       -exp((1/n).*sum(cos(2*pi*x),2)));
 return