一.产生背景

   

❃粒子群算法(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公式

python UAV 粒子群算法 粒子群算法案例_TSP问题

(3)基本PSO算法流程图

python UAV 粒子群算法 粒子群算法案例_PSO_02

关于每个粒子的更新速度和位置的公式如下:

python UAV 粒子群算法 粒子群算法案例_粒子群算法_03

三.简单应用

  

python UAV 粒子群算法 粒子群算法案例_粒子群算法求解TSP问题_04

 



(1)•编码:因为问题的维数为5,所以每个粒子为5维的实数向量。



(2)•初始化范围:根据问题要求,设定为[-30,30]。根据前面的参数分析,我们知道,可以将最大速度设定为Vmax=60。



(3)•种群大小:为了说明方便,这里采用一个较小的种群规模,m=5。



(4)•停止准则:设定为最大迭代次数100次。



(5)•惯性权重:采用固定权重0.5。

(6)邻域拓扑结构:使用星形拓扑结构,即全局版本的粒子群优化算法


 

算法执行的过程如下:

python UAV 粒子群算法 粒子群算法案例_粒子群算法求解TSP问题_05

python UAV 粒子群算法 粒子群算法案例_粒子群算法_06

python UAV 粒子群算法 粒子群算法案例_粒子群算法求解TSP问题_07

python UAV 粒子群算法 粒子群算法案例_粒子群算法求解TSP问题_08

python UAV 粒子群算法 粒子群算法案例_粒子群算法_09

 

 


 


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
 
运行结果如下:


 



 




python UAV 粒子群算法 粒子群算法案例_粒子群算法_10



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
 
	}
}



python UAV 粒子群算法 粒子群算法案例_TSP问题_11



python UAV 粒子群算法 粒子群算法案例_粒子群算法_12



python UAV 粒子群算法 粒子群算法案例_python UAV 粒子群算法_13


 


算法代码地址:(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. 
        
    
    

      }