基于粒子群优化的函数优化
1. 粒子群优化算法
1.1 算法思想
粒子群优化算法(PSO:Particle swarm optimization) 是一种进化计算技术(evolutionary computation)。源于对鸟群捕食的行为研究。粒子群优化算法的基本思想:是通过群体中个体之间的协作和信息共享来寻找最优解。
1.2 基本原理
粒子群算法通过设计一种无质量的粒子来模拟鸟群中的鸟,粒子仅具有两个属性:速度和位置,速度代表移动的快慢,位置代表移动的方向。每个粒子在搜索空间中单独的搜寻最优解,并将其记为当前个体极值,并将个体极值与整个粒子群里的其他粒子共享,找到最优的那个个体极值作为整个粒子群的当前全局最优解,粒子群中的所有粒子根据自己找到的当前个体极值和整个粒子群共享的当前全局最优解来调整自己的速度和位置。
1.3 算法流程
-
初始化
首先,我们设置最大迭代次数,目标函数的自变量个数,粒子的最大速度,位置信息为整个搜索空间,我们在速度区间和搜索空间上随机初始化速度和位置,设置粒子群规模,每个粒子随机初始化一个速度。 -
个体极值和全局最优值
定义适应度函数,个体极值为每个粒子找到的最优解,从这些最优解找到一个全局值,叫做本次全局最优解。与历史全局最优比较,进行更新。 -
粒子速度和位置的更新
粒子速度的更新:
v i d k + 1 = w v i d k + c 1 r a n d ( ) ( p i d − x i d k ) + c 2 r a n d ( ) ( p g b e s t − x i d k ) v^{k+1}_{id}=wv^{k}_{id}+c_1rand()(p_{id}-x^{k}_{id})+c_2rand()(p_{gbest}-x^{k}_{id}) vidk+1=wvidk+c1rand()(pid−xidk)+c2rand()(pgbest−xidk)
粒子位置的更新:
x i d k + 1 = x i d k + v k + 1 i d i = 1 , 2 , . . . , m ; d = 1 , 2 , . . . , D ; x^{k+1}_{id}=x^{k}_{id}+v{k+1}_{id} i=1,2,...,m;d=1,2,...,D; xidk+1=xidk+vk+1idi=1,2,...,m;d=1,2,...,D; -
终止条件
a. 设定迭代次数
b. 代数之间的差值满足最小界限
2. 粒子群优化实现函数优化
2.1 代码实现
%% 清空环境
clc
clear
%% 参数初始化
%粒子群算法中的两个参数
c1 = 1.49445;
c2 = 1.49445;
maxgen=500; % 进化次数
sizepop=100; %种群规模
Vmax=1;
Vmin=-1;
popmax=5;
popmin=-5;
%% 产生初始粒子和速度
for i=1:sizepop
%随机产生一个种群
pop(i,:)=5*rands(1,2); %初始种群
V(i,:)=rands(1,2); %初始化速度
%计算适应度
fitness(i)=fun(pop(i,:)); %染色体的适应度
end
%% 个体极值和群体极值
[bestfitness bestindex]=min(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,:)+0.5*V(j,:);
pop(j,find(pop(j,:)>popmax))=popmax;
pop(j,find(pop(j,:)<popmin))=popmin;
if rand>0.98
pop(j,:)=rands(1,2);
end
%适应度值
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
%% 结果分析
plot(yy)
title('最优个体适应度','fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);
2.2 初始状态下的运行结果
—— | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 平均值 |
---|---|---|---|---|---|---|---|---|---|---|---|
运行时间 | 4.553 s | 4.790 s | 4.660 s | 4.872 s | 4.693 s | 4.629 s | 4.551 s | 4.621 s | 4.775 s | 4.461 s | 4.6605 s |
最优适应度 | 1.989918 | 0.045141 | 4.974795 | 1.989918 | 4.974795 | 9.949586 | 3.979836 | 5.969754 | 2.984877 | 4.974795 | 4.183342 |
3. 粒子群优化算法中的参数讨论
3.1 对于惯性权重w线性变化的讨论
3.1.1 代码实现
下面仅给出有变化的代码段,在原始代码的基础上替换这两部分即可。
%% 参数初始化
%粒子群算法中的三个参数
c1 = 1.49445;%加速因子
c2 = 1.49445;
w = 0.8; %惯性权重
w0 = 0.2;
ww = w-w0;
maxgen=1000; % 进化次s数
sizepop=200; %种群规模
Vmax=1; %限制速度围
Vmin=-1;
popmax=5; %变量取值范围
popmin=-5;
dim=10; %适应度函数维数
func=1; %选择待优化的函数,1为Rastrigin,2为Schaffer,3为Griewank
%% 迭代寻优
for i=1:maxgen
w = w - ww/maxgen;
fprintf('第%d代,',i);
fprintf('W:%f,',w);
fprintf('最优适应度%f\n',fitnessgbest);
for j=1:sizepop
%速度更新
V(j,:) = w*V(j,:) + c1*rand*(pbest(j,:) - pop(j,:)) + c2*rand*(gbest - pop(j,:)); %根据个体最优pbest和群体最优gbest计算下一时刻速度
V(j,find(V(j,:)>Vmax))=Vmax; %限制速度不能太大
V(j,find(V(j,:)<Vmin))=Vmin;
%种群更新
pop(j,:)=pop(j,:)+0.5*V(j,:); %位置更新
pop(j,find(pop(j,:)>popmax))=popmax;%坐标不能超出范围
pop(j,find(pop(j,:)<popmin))=popmin;
end
end
3.1.2 运行结果
—— | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 平均值 |
---|---|---|---|---|---|---|---|---|---|---|---|
运行时间 | 5.054 s | 5.043 s | 4.900 s | 5.164 s | 4.907 s | 5.063 s | 4.866 s | 5.118 s | 5.185 s | 5.668 s | 5.0968 s |
最优适应度 | 2.984879 | 2.984877 | 5.053058 | 0.000000 | 0.000000 | 2.199570 | 0.000850 | 1.002319 | 4.974924 | 5.969754 | 2.5170231 |
- 第一次运行结果
- 第二次运行结果
- 第三次运行结果
- 第四次运行结果
- 第五次运行结果
- 第六次运行结果
- 第七次运行结果
- 第八次运行结果
- 第九次运行结果
- 第十次运行结果
(Tips:之后不再进行冗余的运行结果截图,只进行表格统计及分析)
3.1.3 运行结果分析
惯性权重w线性变化的思想:
惯性权重w是对粒子自身运动状态的信任,当初始迭代时,惯性权重w比较大,反应在速度v的计算公式上可以看出初始迭代的时候粒子的速度比较大,具有很好的全局搜索能力,而局部搜索能力较弱。随着迭代次数的累加,w的值越来越小,粒子的速度也越来越小,此时粒子具有很好的局部搜索能力,而全局搜索能力较弱。但是由于斜率恒定,所以速度的改变总是保持同一水平。如果初始迭代的时候并没有产生较好的点,那么随着迭代的累加以及速度的迅速衰减很可能导致最后陷入局部最优值。
将这十次实验结果求取平均值与初始情况下的参数设置进行对比,可以发现运行时间由4.6605s增加到5.0968s,运行效率降低。思考:因为函数中添加了了每次迭代过程中更新惯性权重w的值,使得算法过程中增加了部分计算,由此运行时间增加,效率略显降低。但是,在效率略降低的情况下,最优适应度得到了理想的值,即0,在算法迭代的后期增加了粒子的局部搜索能力,但仍存在陷入局部最优值的可能。
3.2 对于加速常数c1和c2线性变化的讨论
3.2.1 代码实现
%% 参数初始化
%粒子群算法中的三个参数
c1max = 2;
c1min = 1;
c2max = 2;
c2min = 1;
w = 0.8;
maxgen=1000; % 进化次s数
sizepop=200; %种群规模
Vmax=1; %限制速度围
Vmin=-1;
popmax=5; %变量取值范围
popmin=-5;
dim=10; %适应度函数维数
func=1; %选择待优化的函数,1为Rastrigin,2为Schaffer,3为Griewank
%% 迭代寻优
for i=1:maxgen
c1 = c1max - ((c1max-c1min)/maxgen)*i;
c2 = c2max - ((c2max-c2min)/maxgen)*i;
fprintf('第%d代,',i);
fprintf('c1:%f',c1);
fprintf('c2:%f',c2);
fprintf('最优适应度%f\n',fitnessgbest);
for j=1:sizepop
%速度更新
V(j,:) = w*V(j,:) + c1*rand*(pbest(j,:) - pop(j,:)) + c2*rand*(gbest - pop(j,:)); %根据个体最优pbest和群体最优gbest计算下一时刻速度
V(j,find(V(j,:)>Vmax))=Vmax; %限制速度不能太大
V(j,find(V(j,:)<Vmin))=Vmin;
end
end
3.2.2 运行结果
—— | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 平均值 |
---|---|---|---|---|---|---|---|---|---|---|---|
运行时间 | 5.084 s | 5.272 s | 5.318 s | 4.857 s | 4.963 s | 5.026 s | 5.340 s | 4.947 s | 4.969 | 5.303 s | 5.0587 s |
最优适应度 | 3.979836 | 0.000002 | 3.979836 | 0.998331 | 0.995027 | 8.954632 | 1.989933 | 3.979836 | 5.969754 | 1.024138 | 3.1871325 |
3.2.3 运行结果分析
加速常数cl和c2决定了微粒本身经验信息和其他微粒的经验信息对微粒运行轨迹的影响,反映了微粒群之间的信息交流。设置c1较大的值,会使微粒过多地在局部范围内徘徊,而较大的c2的值,则又会促使微粒过早收敛到局部最小值。
本次参数设置将这十次实验结果求取平均值与初始情况下的参数设置进行对比,可以发现运行时间由4.6605s增加到5.0587s,运行效率降低。思考:因为函数中添加了了每次迭代过程中更新惯性权重w的值,使得算法过程中增加了部分计算,由此运行时间增加,效率略显降低。但是,在效率略降低的情况下,最优适应度得到了较为理想的值,并未得到理想最优值0,但相对于初始参数得到的最优适应度值较为理想。
3.3 对于惯性权重w和加速常数c1、c2同时发生线性变化的讨论
3.3.1 代码实现
%% 参数初始化
%粒子群算法中的三个参数
c1max = 2;
c1min = 1;
c2max = 2;
c2min = 1;
wmax = 0.8; %惯性权重
wmin = 0.2;
maxgen=1000; % 进化次s数
sizepop=200; %种群规模
Vmax=1; %限制速度围
Vmin=-1;
popmax=5; %变量取值范围
popmin=-5;
dim=10; %适应度函数维数
func=1; %选择待优化的函数,1为Rastrigin,2为Schaffer,3为Griewank
%% 迭代寻优
for i=1:maxgen
w = wmax - ((wmax-wmin)/maxgen)*i;
c1 = c1max - ((c1max-c1min)/maxgen)*i;
c2 = c2max - ((c2max-c2min)/maxgen)*i;
fprintf('第%d代,',i);
fprintf('W:%f,',w);
fprintf('c1:%f',c1);
fprintf('c2:%f',c2);
fprintf('最优适应度%f\n',fitnessgbest);
for j=1:sizepop
%速度更新
V(j,:) = w*V(j,:) + c1*rand*(pbest(j,:) - pop(j,:)) + c2*rand*(gbest - pop(j,:)); %根据个体最优pbest和群体最优gbest计算下一时刻速度
V(j,find(V(j,:)>Vmax))=Vmax; %限制速度不能太大
V(j,find(V(j,:)<Vmin))=Vmin;
end
end
3.3.2 运行结果
—— | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 平均值 |
---|---|---|---|---|---|---|---|---|---|---|---|
运行时间 | 4.878 s | 4.937 s | 4.851 s | 4.995 s | 4.999 s | 5.247 s | 5.063 s | 4.978 s | 4.969 s | 5.072 s | 4.9989 |
最优适应度 | 0.000000 | 3.979836 | 0.000000 | 0.000000 | 0.994959 | 0.000000 | 2.984877 | 2.005763 | 0.000000 | 0.000000 | 0.996544 |
3.3.3 运行结果分析
本次实验在每次迭代中,惯性权重w和加速常数c1、c2同时发生线性变化,统计十次的运行结果求取平均值。
相对于单独对惯性权重和加速常数两个参数进行线性变化,发现运行效率略低,而当同时发生线性变化时,运行效率相对有所提高,但对于初始参数来说还是加大了计算量,运行效率还是比初始参数的要低一点,但相差不大。
在本次实验中,多次运行结果中得到最优适应度值为0的结果,虽然得到了理想值,但是在观测中,最优适应值较快地发生收敛,容易陷入局部最优,思考应为两个参数同时发生线性变化时,效果并不会很好。
3.4 维数dim和种群规模sizepop的组合讨论
- dim = 20 sizepop = 200
此次运行多次取其中一次实验结果进行分析。
第1000代,W:0.200000,c1:1.000000c2:2.000000最优适应度5.969754
当把维数由10维增加到20维时,最优适应度较快地发生收敛,得到的最优适应度也并不理想。
可以得到结论:
当维数增加时,会陷入局部最优,得到的结果并不好。 之后的实验将降低维数。
- dim = 5 sizepop = 200
第1000代,W:0.200000,c1:1.000000c2:2.000000最优适应度0.000000
由运行结果可以看到,当维数Dim由10维降低到5维时,算法得到了理想的最有适应度0,且运行时间相比于之前的实验较快,可以得到结论:当维数Dim降低,算法的计算复杂度降低,也可以较快地得到理想最优值。
- dim = 2 sizepop = 200
第1000代,W:0.200000,c1:1.000000c2:2.000000最优适应度0.000000
当维数Dim降低到2维时,其运行效率和得到最优适应度0的迭代次数和5为的相差不大,思考:当维数变化的区间接近10维时,其运行效率和优化求解效率相差较大。
- dim = 10 sizepop = 300
第1000代,W:0.200000,c1:1.000000c2:2.000000最优适应度0.000023
由运行结果可以看到:当维数为10保持不变,更改种群规模sizepop 的值为300时,其运行效率较低,但优化求解得到了0.000023,并不是理想的优化求解值0,并且在在第100次迭代到500次的迭代中陷入了局部最优,发现以200为梯度更改种群规模并不能得到较好的结果,所以加大更改区间进行求解。
- dim = 10 sizepop = 500
第1000代,W:0.200000,c1:1.000000c2:2.000000最优适应度0.000000
当更改种群规模sizepop 的值为500时,其运行效率大幅度降低,甚至超过了10s,但优化求解效率较高,在迭代次数为400的时候得到了理想的优化求解值0,思考:只有以较大梯度更改种群规模的值才可以得到较为理想的优化求解效果,但对于运行效率来说,缺点过于明显,所以增大种群规模这一做法只适用于只期望得到最优适应度并不在意求解效率。
4. 粒子群优化算法的优缺点讨论
4.1 粒子群优化算法的优点
初始化时限制较少,对问题信息几乎没有什么依赖。
算法的记忆能力,存储最优解的能力较强。
算法的原理简单,实现起来比较容易,可以由以上粒子群优化算法基本流程中看到,其核心内容就是一个位置更新和速度更新的公式,相比于遗传算法实现起来较为方便,不会特别复杂。
有着很强的群智能算法的特性,粒子群的个体间之共享信息,有利于搜索的进行。
4.1 粒子群优化算法的缺点
算法的收敛速度过快,导致局部搜索能力不强,直接导致搜索精度不高。
无法保证做到全局搜索,经常在某个局部最优收敛,在本次实验中的惯性权重w和加速常数c1、c2同时发生线性变化时表现的较为明显,多次得到的最优适应度并不理想,并且从某一迭代次数开始最优适应度不再发生变化,纤舞了局部最优。
参数的选择对算法的性能影响较大,很难调整到最适宜的参数。在本次实验中可以明确地感受到,对于参数的设置中,有为数不多的几次得到了最有适应度0,但极不稳定,仅出现了几次,并不会持续出现。