devbox
煮酒与君饮
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

粒子群优化算法（PSO）_粒子群优化算法流程

作者：煮酒与君饮 | 2024-07-03 14:49:05

赞

踩

粒子群优化算法流程

目录

一、算法基本流程

1.1算法概述

1.2 算法流程图

二、代码实现

三、结果分析

一、算法基本流程

1.1算法概述

粒子群优化（PSO, particle swarm optimization）算法是计算智能领域，除了蚁群算法，鱼群算法之外的一种群体智能的优化算法，该算法最早由Kennedy和Eberhart在1995年提出的，该算法源自对鸟类捕食问题的研究。

算法步骤：

 （1）初始化所有的微粒，包括它们的位置和速度，把个体的历史最优pBest设为当前位置，群体的最优个体作为当前的gBest；
 （2）每一代的进化中，计算每个粒子的适应度函数值；
 （3）对每个微粒，将它的适应值和它经历过的最好位置pbest的作比较，如果较好，则将其作为当前的最好位置pBest;
 （4）对每个微粒，将它的适应值和全局所经历最好位置gbest的作比较，如果较好，则重新设置gBest的索引号;
 （5）变化微粒的速度和位置;
 （6）如未达到结束条件(通常为足够好的适应值或达到一个预设最大代数Gmax) ,回到（2），否则输出gBest并结束。
 

  在每一次迭代过程中，粒子通过个体极值和群体极值更新自身的速度和位置，更新公式 如下：

 

 

1.2 算法流程图

 

 

 

二、代码实现

PSO.m

  1.  
  2. %% 清空环境
  3. clc
  4. clear
  5.  
  6. %% 参数初始化
  7. %粒子群算法中的三个参数
  8. c1 = 1.49445;%加速因子
  9. c2 = 1.49445;
  10. w=0.8   %惯性权重
  11.  
  12. maxgen=1000;   % 进化次s数  
  13. sizepop=200;   %种群规模
  14.  
  15. Vmax=1;       %限制速度围
  16. Vmin=-1;     
  17. popmax=5;    %变量取值范围
  18. popmin=-5;
  19. dim=10;       %适应度函数维数
  20.  
  21. func=1;       %选择待优化的函数，1为Rastrigin,2为Schaffer,3为Griewank
  22. Drawfunc(func);%画出待优化的函数，只画出二维情况作为可视化输出
  23.  
  24. %% 产生初始粒子和速度
  25. for i=1:sizepop
  26.     %随机产生一个种群
  27.     pop(i,:)=popmax*rands(1,dim);    %初始种群
  28.     V(i,:)=Vmax*rands(1,dim);             %初始化速度
  29.                                      %计算适应度
  30.     fitness(i)=fun(pop(i,:),func);   %粒子的适应度
  31. end
  32.  
  33. %% 个体极值和群体极值
  34. [bestfitness bestindex]=min(fitness);
  35. gbest=pop(bestindex,:);   %全局最佳
  36. pbest=pop;                %个体最佳
  37. fitnesspbest=fitness;     %个体最佳适应度值
  38. fitnessgbest=bestfitness; %全局最佳适应度值
  39.  
  40. %% 迭代寻优
  41. for i=1:maxgen
  42.     
  43.     fprintf('第%d代，',i);
  44.     fprintf('最优适应度%f\n',fitnessgbest);
  45.     for j=1:sizepop
  46.         
  47.         %速度更新
  48.         V(j,:) = w*V(j,:) + c1*rand*(pbest(j,:) - pop(j,:)) + c2*rand*(gbest - pop(j,:)); %根据个体最优pbest和群体最优gbest计算下一时刻速度
  49.         V(j,find(V(j,:)>Vmax))=Vmax;   %限制速度不能太大
  50.         V(j,find(V(j,:)<Vmin))=Vmin;
  51.         
  52.         %种群更新
  53.         pop(j,:)=pop(j,:)+0.5*V(j,:);       %位置更新
  54.         pop(j,find(pop(j,:)>popmax))=popmax;%坐标不能超出范围
  55.         pop(j,find(pop(j,:)<popmin))=popmin;
  56.         
  57.         if rand>0.98                         %加入变异种子，用于跳出局部最优值
  58.             pop(j,:)=rands(1,dim);
  59.         end
  60.         
  61.         %更新第j个粒子的适应度值
  62.         fitness(j)=fun(pop(j,:),func); 
  63.    
  64.     end
  65.     
  66.     for j=1:sizepop
  67.         
  68.         %个体最优更新
  69.         if fitness(j) < fitnesspbest(j)
  70.             pbest(j,:) = pop(j,:);
  71.             fitnesspbest(j) = fitness(j);
  72.         end
  73.         
  74.         %群体最优更新
  75.         if fitness(j) < fitnessgbest
  76.             gbest = pop(j,:);
  77.             fitnessgbest = fitness(j);
  78.         end
  79.     end 
  80.     yy(i)=fitnessgbest;    
  81.         
  82. end
  83. %% 结果分析
  84. figure;
  85. plot(yy)
  86. title('最优个体适应度','fontsize',12);
  87. xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);
  88.

 Rastrigin.m

图像为：

 

  1. function y = Rastrigin(x)
  2. % Rastrigin函数
  3. % 输入x,给出相应的y值,在x = ( 0 , 0 ,…, 0 )处有全局极小点0.
  4. % 编制人：
  5. % 编制日期：
  6. [row,col] = size(x);
  7. if  row > 1 
  8.     error( ' 输入的参数错误 ' );
  9. end
  10. y =sum(x.^2-10*cos(2*pi*x)+10);
  11. %y =-y;

 Schaffer.m

函数写法：

此函数在（0,...,0）处有最大值1，因此不需要取相反数。

图像为：

 

  1. function y=Schaffer(x)
  2.  
  3. [row,col]=size(x);
  4. if row>1
  5.     error('输入的参数错误');
  6. end
  7. y1=x(1,1);
  8. y2=x(1,2);
  9. temp=y1^2+y2^2;
  10. y=0.5-(sin(sqrt(temp))^2-0.5)/(1+0.001*temp)^2;
  11. %y=-y;

Griewank.m

图像为：

 

  1. function y=Griewank(x)
  2. %Griewan函数
  3. %输入x,给出相应的y值,在x=(0,0,…,0)处有全局极小点0.
  4. %编制人：
  5. %编制日期：
  6. [row,col]=size(x);
  7. if row>1
  8.     error('输入的参数错误');
  9. end
  10. y1=1/4000*sum(x.^2);
  11. y2=1;
  12. for h=1:col
  13.     y2=y2*cos(x(h)/sqrt(h));
  14. end
  15. y=y1-y2+1;
  16. %y=-y;

 

fun.m

函数用于计算粒子适应度值，选择待优化的函数，1为Rastrigin,2为Schaffer,3为Griewank

  1. function y = fun(x,label)
  2.  
  3. %x           input           输入粒子 
  4. %y           output          粒子适应度值 
  5. if label==1
  6.     y=Rastrigin(x);
  7. elseif label==2
  8.     y=Schaffer(x);
  9. else
  10.     y= Griewank(x);
  11. end
  12.

fun1.m 

  1. function y = fun1(x)
  2. %函数用于计算粒子适应度值
  3. %x           input           输入粒子 
  4. %y           output          粒子适应度值 
  5.  
  6. y=-20*exp(-0.2*sqrt((x(1)^2+x(2)^2)/2))-exp((cos(2*pi*x(1))+cos(2*pi*x(2)))/2)+20+exp(1);
  7.  
  8. %y=x(1)^2-10*cos(2*pi*x(1))+10+x(2)^2-10*cos(2*pi*x(2))+10;

 

Drawfunc.m

画出待优化的函数，只画出二维情况作为可视化输出 

  1. function Drawfunc(label)
  2.  
  3. x=-5:0.05:5;%41列的向量
  4. if label==1
  5.     y = x;
  6.     [X,Y] = meshgrid(x,y);
  7.     [row,col] = size(X);
  8.     for  l = 1 :col
  9.          for  h = 1 :row
  10.             z(h,l) = Rastrigin([X(h,l),Y(h,l)]);
  11.         end
  12.     end
  13.     surf(X,Y,z);
  14.     shading interp
  15.     xlabel('x1-axis'),ylabel('x2-axis'),zlabel('f-axis'); 
  16.     title('mesh'); 
  17. end
  18.  
  19. if label==2
  20.     y = x;
  21.     [X,Y] = meshgrid(x,y);
  22.     [row,col] = size(X);
  23.     for  l = 1 :col
  24.          for  h = 1 :row
  25.             z(h,l) = Schaffer([X(h,l),Y(h,l)]);
  26.         end
  27.     end
  28.     surf(X,Y,z);
  29.     shading interp 
  30.     xlabel('x1-axis'),ylabel('x2-axis'),zlabel('f-axis'); 
  31.     title('mesh'); 
  32. end
  33.  
  34. if label==3
  35.     y = x;
  36.     [X,Y] = meshgrid(x,y);
  37.     [row,col] = size(X);
  38.     for  l = 1 :col
  39.          for  h = 1 :row
  40.             z(h,l) = Griewank([X(h,l),Y(h,l)]);
  41.         end
  42.     end
  43.     surf(X,Y,z);
  44.     shading interp 
  45.     xlabel('x1-axis'),ylabel('x2-axis'),zlabel('f-axis'); 
  46.     title('mesh'); 
  47. end
  48.  
  49.  
  50.

PSOMutation.m

  1.  
  2. %% 清空环境
  3. clc
  4. clear
  5.  
  6. %% 参数初始化
  7. %粒子群算法中的两个参数
  8. c1 = 1.49445;
  9. c2 = 1.49445;
  10.  
  11. maxgen=500;   % 进化次数  
  12. sizepop=100;   %种群规模
  13.  
  14. Vmax=1;
  15. Vmin=-1;
  16. popmax=5;
  17. popmin=-5;
  18.  
  19. %% 产生初始粒子和速度
  20. for i=1:sizepop
  21.     %随机产生一个种群
  22.     pop(i,:)=5*rands(1,2);    %初始种群
  23.     V(i,:)=rands(1,2);  %初始化速度
  24.     %计算适应度
  25.     fitness(i)=fun(pop(i,:));   %染色体的适应度
  26. end
  27.  
  28. %% 个体极值和群体极值
  29. [bestfitness bestindex]=min(fitness);
  30. zbest=pop(bestindex,:);   %全局最佳
  31. gbest=pop;    %个体最佳
  32. fitnessgbest=fitness;   %个体最佳适应度值
  33. fitnesszbest=bestfitness;   %全局最佳适应度值
  34.  
  35. %% 迭代寻优
  36. for i=1:maxgen
  37.     
  38.     for j=1:sizepop
  39.         
  40.         %速度更新
  41.         V(j,:) = V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));
  42.         V(j,find(V(j,:)>Vmax))=Vmax;
  43.         V(j,find(V(j,:)<Vmin))=Vmin;
  44.         
  45.         %种群更新
  46.         pop(j,:)=pop(j,:)+0.5*V(j,:);
  47.         pop(j,find(pop(j,:)>popmax))=popmax;
  48.         pop(j,find(pop(j,:)<popmin))=popmin;
  49.         
  50.         if rand>0.98     
  51.             pop(j,:)=rands(1,2);
  52.         end
  53.         
  54.         %适应度值
  55.         fitness(j)=fun(pop(j,:)); 
  56.    
  57.     end
  58.     
  59.     for j=1:sizepop
  60.         
  61.         %个体最优更新
  62.         if fitness(j) < fitnessgbest(j)
  63.             gbest(j,:) = pop(j,:);
  64.             fitnessgbest(j) = fitness(j);
  65.         end
  66.         
  67.         %群体最优更新
  68.         if fitness(j) < fitnesszbest
  69.             zbest = pop(j,:);
  70.             fitnesszbest = fitness(j);
  71.         end
  72.     end 
  73.     yy(i)=fitnesszbest;    
  74.         
  75. end
  76. %% 结果分析
  77. plot(yy)
  78. title('最优个体适应度','fontsize',12);
  79. xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);

 

三、结果分析

POS.m运行结果

最优个体适应度最接近于零时算法达到很好

 

算法中的参数：

 

  1. %粒子群算法中的三个参数
  2. 加速因子
  3. c1 向自身推进的加速权值
  4. c2 向全局推进的加速权值
  5. w  惯性权重
  6.  
  7. sizepop 种群规模
  8.  
  9. dim 适应度函数维数

 

初始时设置：c1 = 1.49445，c2 = 1.49445，w = 0.8，sizepop=200，dim=10

（1）c1 = 1.49445，c2 = 1.49445，sizepop=200，dim=10，改变惯性权重的值  w

w适应度平均值
0.83.581852
0.72.586893
0.62.967745
0.41.893848

w关系到前一速度对当前速度的影响，一般取值在0.9~0.4。因初始种群是随机产生的，当权重w固定为某一个值的时候，我们可以用多次结果取平均值的方法来处理数据。

（2）利用线性递减的函数设置w从0.9到0.4变化，观察记录这时多个结果不同

 典型线性递减策略的w计算公式如下：

w=wmax-(wmax-wmin)\times t\div maxgen

其中 wmax是惯性权重最大值，wmin是惯性权重最小值，t 是当前迭代次数，是个变量，

maxgen是总共迭代次数，是常量，实验中取 1000。

当初始迭代时，惯性权重w比较大，反应在速度v的计算公式上（V(j,:) = w*V(j,:) + c1*rand*(pbest(j,:) - pop(j,:)) + c2*rand*(gbest - pop(j,:))）可以看出初始迭代的时候粒子的速度比较大，具有很好的全局搜索能力，而局部搜索能力较弱。随着迭代次数的累加，w的值越来越小，粒子的速度也越来越小，此时粒子具有很好的局部搜索能力，而全局搜索能力较弱。但是由于斜率恒定，所以速度的改变总是保持同一水平。

 12345678910平均适应度
13.9798363.97983601.9899182.5383814.9747956.96471305.9697544.9747953.537203
200.0000090.9949625.9697543.979836004.9747954.9771370.9950122.189151
30.9949596.9647134.9747951.0481674.9747955.9697540.9949593.9914312.9848772.9848773.588333

 

在迭代寻找最优的代码循环里利用公式来取惯性权重w：

  1. %% 迭代寻优
  2. for i=1:maxgen
  3.     
  4.     fprintf('第%d代，',i);
  5.     fprintf('最优适应度%f\n',fitnessgbest);
  6.     for j=1:sizepop
  7.         
  8.         %速度更新
  9.         V(j,:) = w*V(j,:) + c1*rand*(pbest(j,:) - pop(j,:)) + c2*rand*(gbest - pop(j,:)); %根据个体最优pbest和群体最优gbest计算下一时刻速度
  10.         V(j,find(V(j,:)>Vmax))=Vmax;   %限制速度不能太大
  11.         V(j,find(V(j,:)<Vmin))=Vmin;
  12.         
  13.         %种群更新
  14.         pop(j,:)=pop(j,:)+0.5*V(j,:);       %位置更新
  15.         pop(j,find(pop(j,:)>popmax))=popmax;%坐标不能超出范围
  16.         pop(j,find(pop(j,:)<popmin))=popmin;
  17.         
  18.         if rand>0.98                         %加入变异种子，用于跳出局部最优值
  19.             pop(j,:)=rands(1,dim);
  20.         end
  21.         
  22.         %更新第j个粒子的适应度值
  23.         fitness(j)=fun(pop(j,:),func); 
  24.    
  25.     end
  26.     
  27.     for j=1:sizepop
  28.         
  29.         %个体最优更新
  30.         if fitness(j) < fitnesspbest(j)
  31.             pbest(j,:) = pop(j,:);
  32.             fitnesspbest(j) = fitness(j);
  33.         end
  34.         
  35.         %群体最优更新
  36.         if fitness(j) < fitnessgbest
  37.             gbest = pop(j,:);
  38.             fitnessgbest = fitness(j);
  39.         end
  40.     end 
  41.     yy(i)=fitnessgbest;   
  42.     w=wmax-(wmax-wmin)*i/1000;
  43.         
  44. end

 

 线性微分递减策略的w的计算公式如下：

w=wmax-(wmax-wmin)\times t^{2}\div maxgen^{2}

初始迭代的时候，w变化缓慢，有利于在初始迭代时寻找满足条件的局部最优值，在接近最大迭代次数时，w变化较快，在寻找到局部最优值之后能够快速地收敛逼近于全局最优值，提高运算效率。

（3）其他参数不变，只改变适应度函数维度dim：

dim最优适应度收敛（次）
50469
60697
80752
100483

 

  1. pop(i,:)=popmax*rands(1,dim);    %初始种群
  2. V(i,:)=Vmax*rands(1,dim);             %初始化速度

适应度函数维度用于初始化种群和速度，以及后面若是算法陷入局部最优时加入的变异因子。当dim越大，初始种群和速度比较大的几率高。实验取dim=5或者dim=10时迭代收敛得比较快。 

（4）其他参数不变，只c1、c2改变

适应度c2=1.08324c2=1.49445c2=1.99857
c1=1.083243.9798360.9949593.979836
c1=1.494453.9798361.9899184.974795
c1=1.998575.9697632.9848771.989918

可以看出在此目标函数中，相对三个取值中，c1取1.08324，c2取1.49445，适应度的值更接近于0。

（5）其他参数不变，改变种群规模sizepop

sizepop最优适应度达到最优适应度时收敛次数总用时
10003756.041 s
20005308.334 s
300011512.173 s
40008714.392 s
500052318.376 s
600023920.882 s

   种群规模sizepop影响算法的搜索能力和计算量。当种群规模越大时，达到最优适应度所用的时间越长，此目标函数中种群规模设置在300或400时达到最优适应度时收敛次数比较小，算法效果较好。

 

 

 

 

 

 

 

 

 

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/煮酒与君饮/article/detail/783533
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号