粒子群算法

一.产生背景

   

❃粒子群算法(particleswarm optimization，PSO)由Kennedy和Eberhart在1995年提出，该算法对于Hepper的模拟鸟群(鱼群)的模型进行修正，以使粒子能够飞向解空间，并在最好解处降落，从而得到了粒子群优化算法。

❃同遗传算法类似，也是一种基于群体叠代的，但并没有遗传算法用的交叉以及变异，而是粒子在解空间追随最优的粒子进行搜索。

❃PSO的优势在于简单，容易实现，无需梯度信息，参数少，特别是其天然的实数编码特点特别适合于处理实优化问题。同时又有深刻的智能背景，既适合科学研究，又特别适合工程应用。

设想这样一个场景：一群鸟在随机的搜索食物。在这个区域里只有一块食物，所有的鸟都不知道食物在哪。但是它们知道自己当前的位置距离食物还有多远。

                         那么找到食物的最优策略是什么？

最简单有效的就是搜寻目前离食物最近的鸟的周围区域。

二.算法介绍

(1)简述

❃每个寻优的问题解都被想像成一只鸟，称为“粒子”。所有粒子都在一个D维空间进行搜索。

❃所有的粒子都由一个fitness-function确定适应值以判断目前的位置好坏。

❃每一个粒子必须赋予记忆功能，能记住所搜寻到的最佳位置。

❃每一个粒子还有一个速度以决定飞行的距离和方向。这个速度根据它本身的飞行经验以及同伴的飞行经验进行动态调整。 

(2) 基本PSO算法

  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)

  PS:一般来说，粒子的位置和速度都是在连续的实数空间内进行取值。

   b.基本PSO公式

(3)基本PSO算法流程图

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

三.简单应用

  

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

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

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

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

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

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

算法执行的过程如下:

四.代码实现:运用粒子群算法解决 TSP问题

1.matlab实现

[plain]  view plain  copy

1. close all;  
2. clear all;  
3.   
4. PopSize=500;%种群大小  
5. CityNum = 14;%城市数  
6.   
7. OldBestFitness=0;%旧的最优适应度值  
8.   
9. Iteration=0;%迭代次数  
10. MaxIteration =2000;%最大迭代次数  
11. IsStop=0;%程序停止标志   
12. Num=0;%取得相同适应度值的迭代次数  
13.   
14. c1=0.5;%认知系数  
15. c2=0.7;%社会学习系数  
16. w=0.96-Iteration/MaxIteration;%惯性系数,随迭代次数增加而递减  
17.   
18. %节点坐标  
19. node=[16.47 96.10; 16.47 94.44; 20.09 92.54; 22.39 93.37; 25.23 97.24;...  
20.      22.00 96.05; 20.47 97.02; 17.20 96.29; 16.30 97.38; 14.05 98.12;...  
21.      16.53 97.38; 21.52 95.59; 19.41 97.13; 20.09 94.55];  
22.   
23. %初始化各粒子，即产生路径种群  
24. Group=ones(CityNum,PopSize);     
25. for i=1:PopSize  
26.     Group(:,i)=randperm(CityNum)';  
27. end  
28. Group=Arrange(Group);  
29.   
30. %初始化粒子速度（即交换序）  
31. Velocity =zeros(CityNum,PopSize);     
32. for i=1:PopSize  
33.     Velocity(:,i)=round(rand(1,CityNum)'*CityNum); %round取整  
34. end  
35.   
36. %计算每个城市之间的距离  
37. CityBetweenDistance=zeros(CityNum,CityNum);     
38. for i=1:CityNum  
39.     for j=1:CityNum  
40.         CityBetweenDistance(i,j)=sqrt((node(i,1)-node(j,1))^2+(node(i,2)-node(j,2))^2);  
41.     end  
42. end  
43.   
44. %计算每条路径的距离  
45. for i=1:PopSize     
46.         EachPathDis(i) = PathDistance(Group(:,i)',CityBetweenDistance);  
47. end  
48.   
49. IndivdualBest=Group;%记录各粒子的个体极值点位置,即个体找到的最短路径  
50. IndivdualBestFitness=EachPathDis;%记录最佳适应度值,即个体找到的最短路径的长度  
51. [GlobalBestFitness,index]=min(EachPathDis);%找出全局最优值和相应序号   
52.   
53. %初始随机解  
54. figure;  
55. subplot(2,2,1);  
56. PathPlot(node,CityNum,index,IndivdualBest);  
57. title('随机解');  
58.   
59. %寻优  
60. while(IsStop == 0) & (Iteration < MaxIteration)   
61.     %迭代次数递增  
62.     Iteration = Iteration +1;    
63.       
64.     %更新全局极值点位置,这里指路径  
65.     for i=1:PopSize     
66.         GlobalBest(:,i) = Group(:,index);  
67.         
68.     end  
69.       
70.     %求pij-xij ,pgj-xij交换序，并以概率c1，c2的保留交换序  
71.     pij_xij=GenerateChangeNums(Group,IndivdualBest);    
72.     pij_xij=HoldByOdds(pij_xij,c1);   
73.     pgj_xij=GenerateChangeNums(Group,GlobalBest);  
74.     pgj_xij=HoldByOdds(pgj_xij,c2);  
75.       
76.     %以概率w保留上一代交换序  
77.     Velocity=HoldByOdds(Velocity,w);  
78.   
79.     Group = PathExchange(Group,Velocity); %根据交换序进行路径交换  
80.     Group = PathExchange(Group,pij_xij);  
81.     Group = PathExchange(Group,pgj_xij);  
82.     for i = 1:PopSize    % 更新各路径总距离  
83.           EachPathDis(i) = PathDistance(Group(:,i)',CityBetweenDistance);  
84.       
85.     end  
86.   
87.     IsChange = EachPathDis<IndivdualBestFitness;%更新后的距离优于更新前的，记录序号  
88.     IndivdualBest(:, find(IsChange)) = Group(:, find(IsChange));%更新个体最佳路径  
89.     IndivdualBestFitness = IndivdualBestFitness.*( ~IsChange) + EachPathDis.*IsChange;%更新个体最佳路径距离  
90.     [GlobalBestFitness, index] = min(EachPathDis);%更新全局最佳路径,记录相应的序号  
91.      
92.     if GlobalBestFitness==OldBestFitness %比较更新前和更新后的适应度值;  
93.         Num=Num+1; %相等时记录加一;  
94.     else  
95.         OldBestFitness=GlobalBestFitness;%不相等时更新适应度值，并记录清零;  
96.         Num=0;  
97.     end      
98.     if Num >= 20 %多次迭代的适应度值相近时程序停止  
99.         IsStop=1;  
100.     end  
101.   
102.      BestFitness(Iteration) =GlobalBestFitness;%每一代的最优适应度  
103.   
104.   
105. end  
106.   
107. %最优解  
108. subplot(2,2,2);  
109. PathPlot(node,CityNum,index,IndivdualBest);  
110. title('优化解');  
111. %进化曲线  
112. subplot(2,2,3);  
113. plot((1:Iteration),BestFitness(1:Iteration));  
114. grid on;  
115. title('进化曲线');  
116. %最小路径值  
117. GlobalBestFitness  
118.   
119. 运行结果如下:  

2.java 实现

[java]  view plain  copy

1. package pso;  
2. import java.awt.*;  
3. import java.awt.event.*;  
4. import java.io.ByteArrayInputStream;  
5. import java.io.InputStream;  
6.   
7. import javax.swing.*;  
8. import javax.swing.event.*;  
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 r2;  
21.     private static float x[][];  
22.     private static float v[][];  
23.     private static float xpbest[][];  
24.     private static float pbest[];        
25.     private static float gbest=0;  
26.     private static float xgbest[];  
27.     private static float w;           //惯性因子  
28.     private static float s;  
29.     private static float h;  
30.     private static float fit[];  
31.     public Sounds sound;  
32.       
33.     //粒子群的迭代函数  
34. public void lzqjs()  
35. {  
36.         
37.         w=(float)(0.9-k*(0.9-0.4)/iterations);  
38.         for(int i=0;i<particleNumber;i++)  
39.         {  
40.                    fit[i]= (float)(1/(Math.pow(x[i][0],2)+Math.pow(x[i][1],2))); //求适值函数最大值  
41.                    System.out.print("粒子"+i+"本次适应值函数f为：" + fit[i]);  
42.                    System.out.println();  
43.                    if(fit[i]>pbest[i])  
44.                    {  
45.                     pbest[i]=fit[i];  
46.                     xpbest[i][0]=x[i][0];  
47.                     xpbest[i][1]=x[i][1];  
48.                    }  
49.                    if(pbest[i]>gbest)  
50.                    {  
51.                     gbest=pbest[i];  
52.                     xgbest[0]=xpbest[i][0];  
53.                     xgbest[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.                     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]));  
62.                      
63.                     x[i][j]=(float)(x[i][j]+v[i][j]);  
64.                      
65.                    }  
66.                 System.out.print("粒子"+i+"本次X1的速度变化幅度:"+v[i][0]+";本次X2的速度变化幅度:"+v[i][1]);  
67.                 System.out.println();  
68.                 System.out.print("粒子"+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.         particleNumber=Integer.parseInt(JOptionPane.showInputDialog("请输入粒子个数1-500）"));  
76.         iterations=Integer.parseInt(JOptionPane.showInputDialog("请输入迭代次数"));  
77.         x=new float [particleNumber][2];  
78.         v=new float [particleNumber][2];  
79.         fit=new float [particleNumber];    //存储适值函数值  
80.         pbest=new float [particleNumber];  //存储整个粒子群的最有位置  
81.         xpbest=new float [particleNumber][2];  
82.         xgbest=new float [2];  
83.         for(int i=0;i<particleNumber;i++)  
84.         {  
85.               
86.             //对数组的初始化操作  
87.             pbest[i]=0;  
88.             xpbest[i][0]=0;  
89.             xpbest[i][1]=0;  
90.         }  
91.         xgbest[0]=0;  
92.         xgbest[1]=0;  
93.          System.out.println("开始初始化：");  
94.         for(int i=0;i<particleNumber;i++)  
95.         {  
96.               
97.             for(int j=0;j<2;j++)  
98.             {  
99.                 //任意给定每个位置一定的位置值和速度值  
100.                 x[i][j]=(float)(WMAX*Math.random()+WMIN);  
101.                 v[i][j]=(float)(VMAX*Math.random());  
102.             }  
103.             System.out.print("粒子"+i+"本次X1的变化幅度:"+v[i][0]+";本次X2的变化幅度:"+v[i][1]);  
104.              System.out.println();  
105.             System.out.print("粒子"+i+"本次X1为："+x[i][0]+";本次X2为："+x[i][1]);  
106.              System.out.println();  
107.         }  
108.         System.out.println("初始化数据结束，开始迭代.....");  
109.     Pso threada=new Pso();  
110.     threada.setTitle("基于粒子群的粒子位置动态显示");  
111.     threada.setSize(800,800);  
112.     threada.addWindowListener(new gbck());  
113.     threada.setVisible(true);  
114.         Thread threadc=new Thread(threada);  
115.         threadc.start();  
116.     }  
117.     static class gbck extends WindowAdapter  
118.     {  
119.         public void windowClosing(WindowEvent e)  
120.         {  
121.             System.exit(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.            g.setColor(new Color(0,0,0));  
139.            for(int i=0;i<particleNumber;i++)  
140.            {  
141.             g.drawString("*",(int)(x[i][0]+200),(int)(x[i][1]+200));  
142.            }  
143.            g.setColor(new Color(255,0,0));  
144.            g.drawString("全局最优适应度函数值："+gbest+"      参数1："+xgbest[0]+"     参数2："+xgbest[1]+"    迭代次数："+ k,50,725);  
145.   
146.     try  
147.     {  
148.     lzqjs();  //开始迭代  
149.       
150.     if(k>=iterations)  
151.     {  
152.           
153.         Thread.sleep((int)(5000));  
154.         System.exit(0);  
155.     }  
156.     k=k+1;  //每次迭代一次加1操作  
157.     Thread.sleep((int)(1000));  
158.     }  
159.     catch(InterruptedException e)  
160.     {  
161.          System.out.println(e.toString());  
162.     }  
163.     repaint();  
164.     }  
165.     public  void sound(){  
166.           sound =new Sounds("050.wav");  
167.           InputStream stream =new ByteArrayInputStream(sound.getSamples());  
168.           // play the sound  
169.           sound.play(stream);  
170.           // exit  
171.   
172.     }  
173. }  

运行的结果如下:

算法代码地址:http://download.youkuaiyun.com/detail/u012017783/9700118(Matlab ,java两个版本)