【老生谈算法】matlab编写PSO算法及实例——PSO算法

用MATLAB编写PSO算法及实例

---

1、文档下载：

本算法已经整理成文档如下，有需要的朋友可以点击进行下载

---

2、算法详解：

1.1 粒子群算法
PSO从这种模型中得到启示并用于解决优化问题。PSO 中，每个优化问题的潜在解都是搜索空间中的一只鸟，称之为粒子。所有的粒子都有一个由被优化的函数决定的适值( fitness value) ，每个粒子还有一个速度决定它们飞翔的方向和距离。然后粒子们就追随当前的最优粒子在解空间中搜索。

PSO初始化为一群随机粒子(随机解)，然后通过迭代找到最优解。在每一次迭代中，粒子通过跟踪两个极值来更新自己；第一个就是粒子本身所找到的最优解，这个解称为个体极值；另一个极值是整个种群目前找到的最优解，这个极值是全局极值。另外也可以不用整个种群而只是用其中一部分作为粒子的邻居，那么在所有邻居中的极值就是局部极值。

其中：c1和c2为学习因子，也称加速常数(acceleration constant)，r1和r2为[0，1]范围内的均匀随机数。式(1.1)右边由三部分组成，第一部分为“惯性(inertia)”或“动量(momentum)”部分，反映了粒子的运动“习惯(habit)”，代表粒子有维持自己先前速度的趋势；第二部分为“认知(cognition)”部分，反映了粒子对自身历史经验的记忆(memory)或回忆(remembrance)，代表粒子有向自身历史最佳位置逼近的趋势；第三部分为“社会(social)”部分，反映了粒子间协同合作与知识共享的群体历史经验。

二、算法设计
2.1 算法流程图

2.2 算法实现
算法的流程如下：

2.3 参数选择

（1）种群规模
通常，种群太小则不能提供足够的采样点，以致算法性能很差；种群太大尽管可以增加优化信息，阻止早熟收敛的发生，但无疑会增加计算量，造成收敛时间太长，表现为收敛速度缓慢。种群规模一般设为100~1000。本文选择种群规模为100。

（2）最大迭代次数
迭代次数越多能保证解的收敛性，但是影响运算速度，本文选1000次。

（3）惯性权值
惯性权重表示在多大程度上保留原来的速度。较大，全局收敛能力强，局部收敛能力弱；较小，局部收敛能力强，全局收敛能力弱。本文选0.6。

（4）加速因子
加速常数和分别用于控制粒子指向自身或邻域最佳位置的运动。文献[20]建议，并通常取。本文也取。

（5）粒子维数
本文中粒子维数取决于待优化函数的维数。
需要说明的是，本文的程序允许改变这些参数，因为本文编写的程序参照matlab工具箱，留给用户解决这类问题一个接口函数，上述的各个参数正是接口函数的参数，因此允许改变。另外对于w和c也可采用变参数法，即随迭代次数增加，利用经验公式使它们动态调整，本文采用固定值。

3.1求三维函数f=x(1).^2+x(2).2+x(3).^2 的最小值
步骤：1.初始化x,v; 2.求出每个粒子的适应值；3.初始化pb,pg个体最优和全局最优；4.根据式子更新x,v; 5.是否满足条件，满足跳出循环，否则重复2-4步

尝试编码：
（1）pso.m 文件

%此算法是PSO算法，汪汪的20161024号版本
function[xm,fv]=PSO(fitness,N,c1,c2,w,M,D)
%{
xm,fv算法最后得到的最优解时的x及最优解，fitness为适应度，即要优化的目标函数，
N为种群数量，c1,c2为学习因子，w为惯性权重，M为迭代次数，D为粒子的维数
%}
format long;
%初始化种群
for i=1:N
    for j=1:D
        x(i,j)=randn; %随机初始化位置
        v(i,j)=randn; %随机初始化速度
    end
end
%先计算各个粒子的适应度pi，并初始化y-粒子个体极值，pg-全局极值
for i=1:N
    p(i)=fitness(x(i,:)); %适应度 问题：将x(i,:)改成x(i,j)是否可以，答不能
    y(i,:)=x(i,:); %个体极值
end
pg=x(N,:);%初始化全局极值/最优
for i=1:N-1
    if fitness(x(i,:))<fitness(pg)
        pg=x(i,:);%替换并选出全局极值
    end
end
%进入粒子群算法主要循环，更新v及x
for t=1:M
    for i=1:N
        v(i,:)=w*v(i,:)+c1*rand*(y(1,:)-x(1,:))+c2*rand*(pg-x(i,:));
        x(i,:)=x(i,:)+v(i,:);
        if fitness(x(i,:))<p(i)
            p(i)=fitness(x(i,:));
            y(i,:)=x(i,:);
        end
        if p(i)<fitness(pg)
            pg=y(i,:);
        end
    end
    pbest(t)=fitness(pg);%M次迭代后最优解
end
xm=pg';%为何要共轭转置？
fv=fitness(pg);
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

（2）目标函数fitness.m文件

function f=fitness(x)
f=x(1).^2+x(2).^2+x(3).^2 ;
end
1
2
3

需要说明的是，针对不同的函数优化，只需要改变目标函数就可以。

（3）在命令行输入或建立调用m文件
在命令行先后输入[xm,fv] = PSO(@fitness,100,2,2,0.6,1000,3)，或建立包涵该语句的m文件，运行即可得到结果。

四、 结果与分析
xm=1.0e-04 *
-0.285730140229565
0.676783696397148
-0.250529540096653
fv =
6.024429352056337e-09
fv是最优值，xm为最优值对应的自变量值。

3.2 高斯函数

[x y]=meshgrid(-100:100,-100:100);   %生成xy坐标系且选取范围
sigma=50;    %高斯函数中的sigma为50
img = (1/(2*pi*sigma^2))*exp(-(x.^2+y.^2)/(2*sigma^2)); 
%目标函数，高斯函数
mesh(img);  %三维画图
hold on;  
1
2
3
4
5
6