一.产生背景
❃粒子群算法(particleswarm optimization,PSO)由Kennedy和Eberhart在1995年提出,该算法对于Hepper的模拟鸟群(鱼群)的模型进行修正,以使粒子能够飞向解空间,并在最好解处降落,从而得到了粒子群优化算法。
❃同遗传算法类似,也是一种基于群体叠代的,但并没有遗传算法用的交叉以及变异,而是粒子在解空间追随最优的粒子进行搜索。
❃PSO的优势在于简单,容易实现,无需梯度信息,参数少,特别是其天然的实数编码特点特别适合于处理实优化问题。同时又有深刻的智能背景,既适合科学研究,又特别适合工程应用。
设想这样一个场景:一群鸟在随机的搜索食物。在这个区域里只有一块食物,所有的鸟都不知道食物在哪。但是它们知道自己当前的位置距离食物还有多远。
?
最简单有效的就是搜寻目前离食物最近的鸟的周围区域。
二.算法介绍
(1)简述
❃每个寻优的问题解都被想像成一只鸟,称为“粒子”。所有粒子都在一个D维空间进行搜索。
❃所有的粒子都由一个fitness-function确定适应值以判断目前的位置好坏。
❃每一个粒子必须赋予记忆功能,能记住所搜寻到的最佳位置。
❃每一个粒子还有一个速度以决定飞行的距离和方向。这个速度根据它本身的飞行经验以及同伴的飞行经验进行动态调整。
a. D维空间中,有m个粒子;
粒子i位置:xi=(xi1,xi2,…xiD)
粒子i速度:vi=(vi1,vi2,…viD),1≤i≤m,1 ≤d ≤D
粒子i经历过的历史最好位置:pi=(pi1,pi2,…piD)
群体内(或领域内)所有粒子所经历过的最好位置:
pg =(pg1,pg2,…pgD)
位置和速度都是在连续的实数空间内进行取值。
b.基本PSO公式
(3)基本PSO算法流程图
关于每个粒子的更新速度和位置的公式如下:
三.简单应用
(1)•编码:因为问题的维数为5,所以每个粒子为5维的实数向量。
(2)•初始化范围:根据问题要求,设定为[-30,30]。根据前面的参数分析,我们知道,可以将最大速度设定为Vmax=60。
(3)•种群大小:为了说明方便,这里采用一个较小的种群规模,m=5。
(4)•停止准则:设定为最大迭代次数100次。
(5)•惯性权重:采用固定权重0.5。
(6)邻域拓扑结构:使用星形拓扑结构,即全局版本的粒子群优化算法
算法执行的过程如下:
TSP问题
1.matlab实现
close all;
clear all;
PopSize=500;%种群大小
CityNum = 14;%城市数
OldBestFitness=0;%旧的最优适应度值
Iteration=0;%迭代次数
MaxIteration =2000;%最大迭代次数
IsStop=0;%程序停止标志
Num=0;%取得相同适应度值的迭代次数
c1=0.5;%认知系数
c2=0.7;%社会学习系数
w=0.96-Iteration/MaxIteration;%惯性系数,随迭代次数增加而递减
%节点坐标
node=[16.47 96.10; 16.47 94.44; 20.09 92.54; 22.39 93.37; 25.23 97.24;...
22.00 96.05; 20.47 97.02; 17.20 96.29; 16.30 97.38; 14.05 98.12;...
16.53 97.38; 21.52 95.59; 19.41 97.13; 20.09 94.55];
%初始化各粒子,即产生路径种群
Group=ones(CityNum,PopSize);
for i=1:PopSize
Group(:,i)=randperm(CityNum)';
end
Group=Arrange(Group);
%初始化粒子速度(即交换序)
Velocity =zeros(CityNum,PopSize);
for i=1:PopSize
Velocity(:,i)=round(rand(1,CityNum)'*CityNum); %round取整
end
%计算每个城市之间的距离
CityBetweenDistance=zeros(CityNum,CityNum);
for i=1:CityNum
for j=1:CityNum
CityBetweenDistance(i,j)=sqrt((node(i,1)-node(j,1))^2+(node(i,2)-node(j,2))^2);
end
end
%计算每条路径的距离
for i=1:PopSize
EachPathDis(i) = PathDistance(Group(:,i)',CityBetweenDistance);
end
IndivdualBest=Group;%记录各粒子的个体极值点位置,即个体找到的最短路径
IndivdualBestFitness=EachPathDis;%记录最佳适应度值,即个体找到的最短路径的长度
[GlobalBestFitness,index]=min(EachPathDis);%找出全局最优值和相应序号
%初始随机解
figure;
subplot(2,2,1);
PathPlot(node,CityNum,index,IndivdualBest);
title('随机解');
%寻优
while(IsStop == 0) & (Iteration < MaxIteration)
%迭代次数递增
Iteration = Iteration +1;
%更新全局极值点位置,这里指路径
for i=1:PopSize
GlobalBest(:,i) = Group(:,index);
end
%求pij-xij ,pgj-xij交换序,并以概率c1,c2的保留交换序
pij_xij=GenerateChangeNums(Group,IndivdualBest);
pij_xij=HoldByOdds(pij_xij,c1);
pgj_xij=GenerateChangeNums(Group,GlobalBest);
pgj_xij=HoldByOdds(pgj_xij,c2);
%以概率w保留上一代交换序
Velocity=HoldByOdds(Velocity,w);
Group = PathExchange(Group,Velocity); %根据交换序进行路径交换
Group = PathExchange(Group,pij_xij);
Group = PathExchange(Group,pgj_xij);
for i = 1:PopSize % 更新各路径总距离
EachPathDis(i) = PathDistance(Group(:,i)',CityBetweenDistance);
end
IsChange = EachPathDis<IndivdualBestFitness;%更新后的距离优于更新前的,记录序号
IndivdualBest(:, find(IsChange)) = Group(:, find(IsChange));%更新个体最佳路径
IndivdualBestFitness = IndivdualBestFitness.*( ~IsChange) + EachPathDis.*IsChange;%更新个体最佳路径距离
[GlobalBestFitness, index] = min(EachPathDis);%更新全局最佳路径,记录相应的序号
if GlobalBestFitness==OldBestFitness %比较更新前和更新后的适应度值;
Num=Num+1; %相等时记录加一;
else
OldBestFitness=GlobalBestFitness;%不相等时更新适应度值,并记录清零;
Num=0;
end
if Num >= 20 %多次迭代的适应度值相近时程序停止
IsStop=1;
end
BestFitness(Iteration) =GlobalBestFitness;%每一代的最优适应度
end
%最优解
subplot(2,2,2);
PathPlot(node,CityNum,index,IndivdualBest);
title('优化解');
%进化曲线
subplot(2,2,3);
plot((1:Iteration),BestFitness(1:Iteration));
grid on;
title('进化曲线');
%最小路径值
GlobalBestFitness
运行结果如下:
package pso;
import java.awt.*;
import java.awt.event.*;
import java.io.ByteArrayInputStream;
import java.io.InputStream;
import javax.swing.*;
import javax.swing.event.*;
public class Pso extends Frame implements Runnable
{
private static int particleNumber; //粒子的数量
private static int iterations; //迭代的次数
private static int k=1; //记录迭代的次数
final private static float C1=2; //学习因子
final private static float C2=2;
final private static float WMIN=-200;
final private static float WMAX=200;
final private static float VMAX=200;
private static float r1; //随机数0-1之间
private static float r2;
private static float x[][];
private static float v[][];
private static float xpbest[][];
private static float pbest[];
private static float gbest=0;
private static float xgbest[];
private static float w; //惯性因子
private static float s;
private static float h;
private static float fit[];
public Sounds sound;
//粒子群的迭代函数
public void lzqjs()
{
w=(float)(0.9-k*(0.9-0.4)/iterations);
for(int i=0;i<particleNumber;i++)
{
fit[i]= (float)(1/(Math.pow(x[i][0],2)+Math.pow(x[i][1],2))); //求适值函数最大值
System.out.print("粒子"+i+"本次适应值函数f为:" + fit[i]);
System.out.println();
if(fit[i]>pbest[i])
{
pbest[i]=fit[i];
xpbest[i][0]=x[i][0];
xpbest[i][1]=x[i][1];
}
if(pbest[i]>gbest)
{
gbest=pbest[i];
xgbest[0]=xpbest[i][0];
xgbest[1]=xpbest[i][1];
}
}
for(int i=0;i<particleNumber;i++)
{
for(int j=0;j<2;j++)
{
//粒子速度和位置迭代方程:
v[i][j]=(float)(w*v[i][j]+C1*Math.random()*(xpbest[i][j]-x[i][j])+C2*Math.random()*(xgbest[j]-x[i][j]));
x[i][j]=(float)(x[i][j]+v[i][j]);
}
System.out.print("粒子"+i+"本次X1的速度变化幅度:"+v[i][0]+";本次X2的速度变化幅度:"+v[i][1]);
System.out.println();
System.out.print("粒子"+i+"本次X1为:"+x[i][0]+";本次X2为:"+x[i][1]);
System.out.println();
}
}
public static void main(String[] args)
{
particleNumber=Integer.parseInt(JOptionPane.showInputDialog("请输入粒子个数1-500)"));
iterations=Integer.parseInt(JOptionPane.showInputDialog("请输入迭代次数"));
x=new float [particleNumber][2];
v=new float [particleNumber][2];
fit=new float [particleNumber]; //存储适值函数值
pbest=new float [particleNumber]; //存储整个粒子群的最有位置
xpbest=new float [particleNumber][2];
xgbest=new float [2];
for(int i=0;i<particleNumber;i++)
{
//对数组的初始化操作
pbest[i]=0;
xpbest[i][0]=0;
xpbest[i][1]=0;
}
xgbest[0]=0;
xgbest[1]=0;
System.out.println("开始初始化:");
for(int i=0;i<particleNumber;i++)
{
for(int j=0;j<2;j++)
{
//任意给定每个位置一定的位置值和速度值
x[i][j]=(float)(WMAX*Math.random()+WMIN);
v[i][j]=(float)(VMAX*Math.random());
}
System.out.print("粒子"+i+"本次X1的变化幅度:"+v[i][0]+";本次X2的变化幅度:"+v[i][1]);
System.out.println();
System.out.print("粒子"+i+"本次X1为:"+x[i][0]+";本次X2为:"+x[i][1]);
System.out.println();
}
System.out.println("初始化数据结束,开始迭代.....");
Pso threada=new Pso();
threada.setTitle("基于粒子群的粒子位置动态显示");
threada.setSize(800,800);
threada.addWindowListener(new gbck());
threada.setVisible(true);
Thread threadc=new Thread(threada);
threadc.start();
}
static class gbck extends WindowAdapter
{
public void windowClosing(WindowEvent e)
{
System.exit(0);
}
}
//开启的额外线程用于声音的播放
public void run()
{
repaint();
for(int i=0;i<iterations;i++){
sound();
}
}
public void paint(Graphics g)
{
g.setColor(new Color(0,0,0));
for(int i=0;i<particleNumber;i++)
{
g.drawString("*",(int)(x[i][0]+200),(int)(x[i][1]+200));
}
g.setColor(new Color(255,0,0));
g.drawString("全局最优适应度函数值:"+gbest+" 参数1:"+xgbest[0]+" 参数2:"+xgbest[1]+" 迭代次数:"+ k,50,725);
try
{
lzqjs(); //开始迭代
if(k>=iterations)
{
Thread.sleep((int)(5000));
System.exit(0);
}
k=k+1; //每次迭代一次加1操作
Thread.sleep((int)(1000));
}
catch(InterruptedException e)
{
System.out.println(e.toString());
}
repaint();
}
public void sound(){
sound =new Sounds("050.wav");
InputStream stream =new ByteArrayInputStream(sound.getSamples());
// play the sound
sound.play(stream);
// exit
}
}
算法代码地址:(Matlab ,java两个版本)
1. package
2.
import
3.
import
4.
import
5.
import
6.
7.
import
8.
import
9.
public class Pso extends Frame
implements
Runnable
10.
{
11.
private static int particleNumber; //粒子的数量
12.
private static int iterations; //迭代的次数
13.
private static int k= 1; //记录迭代的次数
14.
final private static float C1= 2;
//学习因子
15.
final private static float C2= 2;
16.
final private static float WMIN=- 200;
17.
final private static float WMAX= 200;
18.
final private static float VMAX= 200;
19.
private static float r1; //随机数0-1之间
20.
private static float
21.
private static float
22.
private static float
23.
private static float
24.
private static float
25.
private static float gbest= 0;
26.
private static float
27.
private static float w; //惯性因子
28.
private static float
29.
private static float
30.
private static float
31.
public
32.
33.
//粒子群的迭代函数
34.
public void lzqjs()
35.
{
36.
37.
float)( 0.9-k*( 0.9- 0.4)/iterations);
38.
for( int i= 0;i<particleNumber;i++)
39.
{
40.
float)( 1/(Math.pow(x[i][ 0], 2)+Math.pow(x[i][ 1],
2)));
//求适值函数最大值
41.
"粒子"+i+ "本次适应值函数f为:"
42.
System.out.println();
43.
if(fit[i]>pbest[i])
44.
{
45.
pbest[i]=fit[i];
46.
0]=x[i][ 0];
47.
1]=x[i][ 1];
48.
}
49.
if(pbest[i]>gbest)
50.
{
51.
gbest=pbest[i];
52.
0]=xpbest[i][ 0];
53.
1]=xpbest[i][ 1];
54.
}
55.
}
56.
for( int i= 0;i<particleNumber;i++)
57.
{
58.
for( int j= 0;j< 2;j++)
59.
{
60.
//粒子速度和位置迭代方程:
61.
float)(w*v[i][j]+C1*Math.random()*(xpbest[i][j]-x[i][j])+C2*Math.random()*(xgbest[j]-x[i][j]));
62.
63.
float)(x[i][j]+v[i][j]);
64.
65.
}
66.
"粒子"+i+ "本次X1的速度变化幅度:"+v[i][ 0]+ ";本次X2的速度变化幅度:"+v[i][ 1]);
67.
System.out.println();
68.
"粒子"+i+ "本次X1为:"+x[i][ 0]+ ";本次X2为:"+x[i][ 1]);
69.
System.out.println();
70.
}
71.
}
72.
public static void main(String[] args)
73.
{
74.
75.
"请输入粒子个数1-500)"));
76.
"请输入迭代次数"));
77.
new float [particleNumber][ 2];
78.
new float [particleNumber][ 2];
79.
new float [particleNumber]; //存储适值函数值
80.
new float [particleNumber]; //存储整个粒子群的最有位置
81.
new float [particleNumber][ 2];
82.
new float [ 2];
83.
for( int i= 0;i<particleNumber;i++)
84.
{
85.
86.
//对数组的初始化操作
87.
0;
88.
0]= 0;
89.
1]= 0;
90.
}
91.
0]= 0;
92.
1]= 0;
93.
"开始初始化:");
94.
for( int i= 0;i<particleNumber;i++)
95.
{
96.
97.
for( int j= 0;j< 2;j++)
98.
{
99.
//任意给定每个位置一定的位置值和速度值
100.
float)(WMAX*Math.random()+WMIN);
101.
float)(VMAX*Math.random());
102.
}
103.
"粒子"+i+ "本次X1的变化幅度:"+v[i][ 0]+ ";本次X2的变化幅度:"+v[i][ 1]);
104.
System.out.println();
105.
"粒子"+i+ "本次X1为:"+x[i][ 0]+ ";本次X2为:"+x[i][ 1]);
106.
System.out.println();
107.
}
108.
"初始化数据结束,开始迭代.....");
109.
new
110.
"基于粒子群的粒子位置动态显示");
111.
800, 800);
112.
new
113.
true);
114.
new
115.
threadc.start();
116.
}
117.
static class gbck extends WindowAdapter
118.
{
119.
public void windowClosing(WindowEvent e)
120.
{
121.
0);
122.
}
123.
}
124.
125.
//开启的额外线程用于声音的播放
126.
public void run()
127.
{
128.
129.
repaint();
130.
131.
for( int i= 0;i<iterations;i++){
132.
sound();
133.
}
134.
}
135.
public void paint(Graphics g)
136.
{
137.
138.
new Color( 0, 0, 0));
139.
for( int i= 0;i<particleNumber;i++)
140.
{
141.
"*",( int)(x[i][ 0]+ 200),( int)(x[i][
1]+
200));
142.
}
143.
new Color( 255, 0, 0));
144.
"全局最优适应度函数值:"+gbest+ " 参数1:"+xgbest[ 0]+ " 参数2:"+xgbest[ 1]+
" 迭代次数:"+ k,
50,
725);
145.
146.
try
147.
{
148.
lzqjs();
//开始迭代
149.
150.
if(k>=iterations)
151.
{
152.
153.
int)( 5000));
154.
0);
155.
}
156.
1; //每次迭代一次加1操作
157.
int)( 1000));
158.
}
159.
catch(InterruptedException e)
160.
{
161.
System.out.println(e.toString());
162.
}
163.
repaint();
164.
}
165.
public void sound(){
166.
new Sounds( "050.wav");
167.
new
168.
// play the sound
169.
sound.play(stream);
170.
// exit
171.
172.
}
173.
}