算法篇----简易版的粒子群（PSO）实现（Matlab语言）

前言

    粒子群算法实现起来并不是很难，算法思想可以参加我上一篇博文，不多说了。好了，Matlab版的粒子群走起。

1 定义变量

    粒子群算法有很多参数，做实验的时候会纠结在参数问题上，这里就随机设定了。有时候，参数好坏，是成败关键。没有修过Matlab语言的朋友不用担心，之后会把C或Java版的也总结一下。知识嘛，多多益善。Let's  go ！

    根据公式，需要c1 c2 r1 r2 w M 

    代码如下：

c1 = 2;%学习因子 c2 = 2;%学习因子 w = 0.7;%惯性权重 D = 10;%搜索维度 M = 40;%种群大小 %r1、r2 由randn生成，不声明了

2 初始化种群

    这里是随机生成在不同维度上的粒子，就是说在D个方向上撒上M个粒子。

    代码如下：

for i=1:M for j=1:D x(i,j)=randn; %随机初始化位置 v(i,j)=randn; %随机初始化速度 end end

3 计算粒子适应度

    通过调用适应度函数，计算适应度值。在此处，初始化个体的适应度值p(i)和全局最优的适应度值gbest。

    代码如下：

for i=1:M p(i)=fitness(x(i,:),D); y(i,:)=x(i,:); end gbest=x(1,:); %gbest为全局最优适应度值

4 更新个体最优值

    将适应度函数计算的结果和目前的个体适应度值相比较，将适应度值较小的位置保存下来。

    代码如下：
for i = 1:M if fitness(x(i,:),D) < p(i) p(i)=fitness(x(i,:),D); y(i,:)=x(i,:); end end

5 更新全局最优值

    上面已经将全局最优位置初始化，只需要将它和目前个体的最优位置相比较即可，位置代入适应度函数，取最小的值对应的位置，保存该位置即可。

    代码如下：

if p(i) < fitness(pg,D) gbest=y(i,:); end

6 更新粒子的速度和位置

    根据公式来写即可，对于不同的粒子群算法这个地方会需要改动。

    代码如下：
for i=1:M v(i,:)=w*v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(pg-x(i,:)); x(i,:)=x(i,:)+v(i,:); end

7 适应度函数

    适应度函数最好另外建一个Script文件，把它放进去。适应度函数可以根据要研究的问题做相应的改动，下面列出一个函数，仅供测试。
    代码如下: function result = fitness(X,D) sum = 0.0; for i = 1:D sum = sum + X(i)^2; end result = sum;

8  完整代码

    到目前，应该有两个m文件，一个是适应度函数的，一个是粒子群的。不过，还需要在粒子群的文件中加入一些代码，主要是调整代码顺序，迭代次数和精度问题，具体参见代码即可。

    完整代码如下：
 clear all; clc; format long; %------给定初始化条件---------------------------------------------- c1=2; %学习因子1 c2=2; %学习因子2 w=0.7; %惯性权重 MaxDT=1000; %最大迭代次数 D=10; %搜索空间维数（未知数个数） M=40; %初始化群体个体数目 eps=10^(-6); %设置精度(在已知最小值时候用) %------初始化种群的个体(可以在这里限定位置和速度的范围)------------ for i=1:M for j=1:D x(i,j)=randn; %随机初始化位置 v(i,j)=randn; %随机初始化速度 end end %------先计算各个粒子的适应度，并初始化p(i)和gbest-------------------- for i=1:M p(i)=fitness(x(i,:),D); y(i,:)=x(i,:); end gbest=x(1,:); %gbest为全局最优 for i=2:M if fitness(x(i,M:),D) < fitness(gbest,D) gbest=x(i,:); end end %------进入主要循环，按照公式依次迭代，直到满足精度要求------------ for t=1:MaxDT for i=1:M v(i,:)=w*v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(gbest-x(i,:)); x(i,:)=x(i,:)+v(i,:); if fitness(x(i,:),D)<p(i) p(i)=fitness(x(i,:),D); y(i,:)=x(i,:); end if p(i)<fitness(gbest,D) gbest=y(i,:); end end end %------显示计算结果 disp('*************************************************************') disp('函数的全局最优位置为：') Solution=gbest' disp('最后得到的优化极值为：') Result=fitness(gbest,D) disp('*************************************************************')

9 实验结果

如图：