算法篇----简易版的粒子群（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 实验结果

如图：