粒子群算法(PSO)详解

原文

一.产生背景

   

粒子群算法(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实现
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
 
运行结果如下:
 
 
2.java 实现
  1. 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
     
    	}
    }
运行的结果如下:

 
算法代码地址:http://download.csdn.net/detail/u012017783/9700118(Matlab ,java两个版本)
  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.  
    }
  • 2
    点赞
  • 8
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值