南山小马

南山小马

0%

【计算智能】目标优化智能算法之粒子群算法

发表于 更新于 分类于 Disqus：

粒子群优化(Particle Swarm Optimization, PSO)，又称微粒算法，是计算智能领域，除了蚁群算法、鱼群算法之外的一种群体智能的优化算法。

历史背景

该算法最早是由J. Kennedy和R. C. Eberhart等于1995年开发的一种演化计算技术。PSO算法源于对鸟类捕食行为的研究，鸟类捕食时，找到食物最简单有效的策略就是搜寻当前距离食物最近的鸟的周围区域。PSO算法是从这种生物种群行为特征中得到启发并用于求解优化问题的，算法中每个粒子都代表问题的一个潜在解，每个粒子对应一个由适应度函数决定的适应度值。粒子的速度决定了粒子的移动方向和距离，速度随自身及其他粒子的移动经验进行动态调整，从而实现个体在可解空间中的寻优。通过加入近邻的速度匹配、并考虑多维搜索和根据距离的加速、形成了PSO的最初版本。之后引入了惯性权重w来更好的控制开发(exploitation)和探索(exploration)，形成了标准版本。为了提高粒子群算的性能和实用性，中山大学、（英国）格拉斯哥大学等又开发了自适应(Adaptive PSO)版本和离散(discrete)版本。

常见应用

在本文的开篇已经提到，粒子群算法是属于智能算法中的目标优化算法，因此，只要与优化问题相关都可以采用粒子群算法。而问题又可以分为连续和离散的，因此粒子群算法也就有连续版本和离散版本。标准粒子群算法是连续的，若要求解的问题是离散的，那么我们就需要根据问题的实际情况，对标准粒子群算法进行修改，以适应离散的情况。
可以参考这篇粒子群算法的综述《粒子群优化算法综述》戴朝华

典型理论问题包括：组合优化、约束优化、多目标优化、动态系统优化等。实际工业应用有：电力系统、滤波器设计、自动控制、数据聚类、模式识别与图像处理、化工、机械、通信、机器人、经济、生物信息、医学、任务分配、TSP等等。

主要思想

PSO算法首先在可行解空间中初始化一群粒子，每个粒子都代表极值优化问题的一个潜在最优解，用位置、速度和适应度值三项指标表示该粒子特征，适应度值由适应度函数计算得到，其值的好坏表示粒子的优劣。粒子在解空间中运动，通过跟踪个体极值pBest和群体极值gBest更新个体位置。个体极值pBest是指个体所经历位置中计算得到的适应度值最优位置，群体极值gBest是指种群中的所有粒子搜索到的适应度最优位置。粒子每更新一次位置，就计算一次适应度值，并且通过比较新粒子的适应度值和个体极值、群体极值的适应度值更新个体极值pBest和群体极值gBest位置。
由于粒子群算法是源于鸟类的捕食行为，那么我们可以想象这么一个场景,以便对粒子群算法有个形象的了解:
在一片广阔的土地上（解空间），有一群鸟（粒子）在觅食，而远程有一片玉米地（目标函数），每一只鸟在飞行（搜索）的过程中，它都记下它吃过最好吃的位置（个体极值），同时它又和其他伙伴交流哪里有好吃的（群体极值），根据这些条件判断，它判断下一时刻应该飞往哪里（循环迭代）。

算法步骤

流程图

  1. 粒子和速度初始化
    如果为了简单起见，可以采用随机的方式进行初始化。但是为了加快算法效率，可以根据问题实际来进行初始化操作。
    假设在一个D维的搜索空间中，由n个粒子组成的种群 $X=({X_2},{X_2},{X_3},\cdots,{X_n})$ 其中第i个粒子表示为一个D为的空间向量 ${X_i}={({x_{i1}},{x_{i2}},\cdots,{x_{id}})^T}$ 代表第i个粒子在D维搜索空间中的位置，亦代表问题的一个潜在解。
  2. 粒子适应度计算
    根据目标函数即可计算出每个粒子位置Xi对应的适应度值。第i个粒子的速度为 ${V_i} = {({V_{i1}},{V_{i2}}, \cdots ,{V_{id}})^T}$ 其个体极值为 ${P_i} = {({P_{i1}},{P_{i2}}, \cdots ,{P_{id}})^T}$ 种群的群体极值为 ${P_g} = {({P_{g1}},{P_{g2}}, \cdots ,{P_{gd}})^T}$
  3. 速度更新和位置更新
    在每次迭代过程中，粒子通过个体极值和群体极值更新自身的速度和位置，即 $$V_{id}^{k + 1} = \omega V_{id}^k + {c_1}{r_1}(P_{id}^k - X_{id}^k) + {c_2}{r_2}(P_{gd}^k - X_{id}^k)$$ $$X_{id}^{k + 1} = X_{id}^k + {V_{k + 1{\kern 1pt} {\kern 1pt} {\kern 1pt} id}}$$ 其中，ω为惯性权重；d=1,2,…,D; i=1,2,…,n; k为当前迭代次数；Vid为粒子的速度；c1和c2是非负的常数，称为加速度因子；r1和r2是分布于[0,1]区间的随机数。为防止粒子的盲目搜索，一遍建议将其位置和速度限制在一定的区间[Xmin,Xmax]、[Vmin,Vmax]

参数设定

  1. 速度更新参数c1和c2 一般设置为1.49445
  2. 惯性权重
    惯性权重表示为粒子本身的飞行速度。顾名思义ω实际反映了粒子过去的运动状态对当前行为的影响，就像是我们物理中提到的惯性。如果ω << 1，从前的运动状态很少能影响当前的行为，粒子的速度会很快的改变；相反，ω较大，虽然会有很大的搜索空间，但是粒子很难改变其运动方向，很难向较优位置收敛，由于算法速度的因素，在实际运用中很少这样设置。也就是说，较高的ω设置促进全局搜索，较低的ω设置促进快速的局部搜索。
    Shi.Y提出了线性递减惯性权重（linear decreasing inertia weight,LDIW），并分析指出一个较大的惯性权值有利于全局搜索，而一个较小的惯性权值则更利于局部搜索。为了更好地平衡算法的全局搜索与局部搜索能力。
    $$\omega (k) = {\omega _{end}} + ({\omega _{start}} - {\omega _{end}})({T_{\max }} - k)/{T_{\max }}$$ 其中，ωstart为初始惯性权重，ωend为迭代到最大次数时的惯性权重；k为当前迭代次数；Tmax为最大迭代次数。一般来说，惯性权值ωstart=0.9，ωend=0.4时算法性能最好。这样，随着迭代的进行，惯性权重有0.9线性递减至0.4，迭代初期较大的惯性权重使算法保持了较强的全局搜索能力，而迭代后期较小的惯性权重有利于算法进行更精确的局部搜索。

MATLAB参考代码

来源《MATLAB智能算法30个案例分析》第十三章，粒子群算法的寻优算法
主函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
%% 清空环境
clc
clear

%% 参数初始化
%粒子群算法中的两个参数
c1 = 1.49445;
c2 = 1.49445;

maxgen=300;    %进化次数  
sizepop=20;   %种群规模
Vmax=0.5;
Vmin=-0.5;
popmax=2;
popmin=-2;

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

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

    %% 迭代寻优
    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
    s(k,:)=yy;
end
%% 结果分析
for m=1:300
    s(101,m)=sum( s(:,m) )/100;
end
plot(s(101,:),'k')
title('最优个体适应度','fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);

适应度函数

1
2
3
4
5
function y = fun(x)
%函数用于计算粒子适应度值
%x           input           输入粒子 
%y           output          粒子适应度值 
y=sin(sqrt(x(1).^2+x(2).^2))./sqrt(x(1).^2+x(2).^2)+exp((cos(2*pi*x(1))+cos(2*pi*x(2)))/2)-2.71289;

参考资料

《MATLAB智能算法30个案例分析》
粒子群优化维基百科
[Algorithm] 群体智能优化算法之粒子群优化算法