粒子群算法

粒子群算法(PSO),在PSO中,每个优化问题的潜在解都是搜索空间的一只鸟，被称为粒子，所有的粒子都有一个由适应度函数决定的适值，每个粒子还有一个速度决定它们“”飞行“”的方向和距离，然后粒子们就追随当前的最优粒子在解空间中搜索,整个过程大致为,PSO初始化为一群随机粒子(随机解),然后通过迭代找到最优解。在每一次迭代的过程中，粒子通过跟踪两个极值来更新自己:第一个就是粒子本身所找到的最优解,这个解称为个体极值;另一个极值是整个种群目前找到的最优解,这个解是全局极值
假设初始目标搜索空间（初始种群位置） D 维 N 列 X i = ( x i 1 , . . . x i D ) i = 1 , . . . N 第 i 个粒子的飞行速度记为 V i = ( v i 1 , . . . v i D ) i = 1 , . . . N 第 i 个粒子的最优位置个体最优记为 P b e s t = ( p i 1 , . . . p i D ) i = 1 , . . . N 全局最优记为 g b e s t = ( P g 1 , . . . p g D ) 根据以下公式更新位置和速度 V i d = W ∗ V i d + c 1 r 1 ( P i d − X i d ) + c 2 r 2 ( P g d − X i d ) X i d = x i d + v i d 其中 w 为惯性银子， c 1 和 c 2 为学习因子， r 1 和 r 2 为随机数 v i d 为粒子当前速度， P i d 为该粒子当前最优位置， P g d 为全局最优位置 \huge假设初始目标搜索空间（初始种群位置）D维N列\X_i=(x_{i1},…x_{iD})i=1,…N\\huge第i个粒子的飞行速度记为\V_i=(v_{i1},…v_{iD})i=1,…N\\huge第i个粒子的最优位置 个体最优记为\P_{best}=(p_{i1},…p_{iD})i=1,…N\\huge全局最优记为\g_{best}=(P_{g1},…p_{gD}) \ \ \ \huge根据以下公式更新位置和速度\V_{id}=W/*V_{id}+c_1r_1(P_{id}-X_{id})+c_2r_2(P_{gd}-X_{id})\X_{id}=x_{id}+v_{id}\\huge其中w为惯性银子，c_1和c_2为学习因子，r_1和r_2为随机数\\huge v_{id}为粒子当前速度，\\huge P_{id}为该粒子当前最优位置，P_{gd}为全局最优位置假设初始目标搜索空间（初始种群位置）D维N列Xi=(xi1,…xiD)i=1,…N第i个粒子的飞行速度记为Vi=(vi1,…viD)i=1,…N第i个粒子的最优位置个体最优记为Pbest=(pi1,…piD)i=1,…N全局最优记为gbest=(Pg1,…pgD)根据以下公式更新位置和速度Vid=W∗Vid+c1r1(Pid−Xid)+c2r2(Pgd−Xid)Xid=xid+vid其中w为惯性银子，c1和c2为学习因子，r1和r2为随机数vid为粒子当前速度，Pid为该粒子当前最优位置，Pgd为全局最优位置

代码

初始化

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
%% 初始化
clear
clc
f=@(x)x.*sin(x).*cos(2*x)-2*x.*sin(3*x)+3*x.*sin(4*x);
N=20;           %初始化种群个数
d=1;            %可行解维数
ger=120;        %最大迭代次数
limit=[-10,10];   %设置位置参数限制
w=0.8;          %惯性权重
c1=0.5;         %自我学习因子
c2=0.5;         %群体学习因子
figure(1);
ezplot(f,[limit(1),0.01,limit(2)]);

x=limit(1)+(limit(2)-limit(1)).*rand(N,d);  %初始种群的位置
v=rand(N,d);                                %初始种群的速度

初始化最优位置

1
2
3
4
5
6
7
8
9
10
11
%% 计算各个粒子的适应度 初始化p（i）和pg
for i=1:N
  p(i)=f(x(i,:)); %各个粒子最优适应度,因为第一代,个体最优适应度就是其本身的适应度
  y(i,:)=x(i,:); %各个粒子的个体最优位置,因为为第一代,个体最位置就是其本身
end
pg=x(N,:); %随意初始化,pg为全局最优位置
for i=1:N-1 %寻找全局最优位置
  if f(x(i,:))>f(pg)
    pg=x(i,:);
  end
end

关键算法

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
%% obj：求解f(x)=xsinxcos2x-2xsin3x+3xsin4x在[0,50]的最小值
%% 初始化
clear
clc
f=@(x)x.*sin(x).*cos(2*x)-2*x.*sin(3*x)+3*x.*sin(4*x);
N=20;           %初始化种群个数
d=1;            %可行解维数
ger=120;        %最大迭代次数
limit=[-10,10];   %设置位置参数限制
w=0.8;          %惯性权重
c1=0.5;         %自我学习因子
c2=0.5;         %群体学习因子
figure(1);
ezplot(f,[limit(1),0.01,limit(2)]);

x=limit(1)+(limit(2)-limit(1)).*rand(N,d);  %初始种群的位置
v=rand(N,d);                                %初始种群的速度

%% 计算各个粒子的适应度 初始化p（i）和pg
for i=1:N
  p(i)=f(x(i,:)); %各个粒子最优适应度,因为第一代,个体最优适应度就是其本身的适应度
  y(i,:)=x(i,:); %各个粒子的个体最优位置,因为为第一代,个体最位置就是其本身
end
pg=x(N,:); %随意初始化,pg为全局最优位置
for i=1:N-1 %寻找全局最优位置
  if f(x(i,:))>f(pg)
    pg=x(i,:);
  end
end

%% 群体更新
iter=1;

while iter<=ger
    for i=1:N
        v(i,:)=w*v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(pg-x(i,:)); %优化粒子速度
        x(i,:)=x(i,:)+v(i,:); %优化粒子位置
        %防止越界,若越界,则取边界值，没有约束的话去掉这两个if语句
        if x(i,:)<limit(1)
          x(i,:)=limit(1);
        end
        if x(i,:)>limit(2)
          x(i,:)=limit(2)
        end
        if f(x(i,:))<p(i) %更新个体最优
          p(i)=f(x(i,:));
          y(i,:)=x(i,:);
        end
        if p(i)<f(pg) %更新全局最优
          pg=y(i,:);
        end
    end
    iter=iter+1;  
end
hold on
plot(pg,f(pg),'ro')
hold off

作者声明

1
如有问题，欢迎指正！