粒子群算法matlab代码吐血推荐.docx

粒子群算法(1)粒子群算法简介
二、粒子群算法的具体表述
　　 上面罗嗦了半天，那些都是科研工作者写论文的语气，不过，的历史就像上面说的那样。下面通俗的解释算法。
　　 算法就是模拟一群鸟寻找食物的过程，每个鸟就是.中的粒子，也就是我们需要求解问题的可能解，这些鸟在寻找食物的过程中，不停改变自己在空中飞行的位置与速度。大家也可以观察一下，鸟群在寻找食物的过程中，开始鸟群比较分散，逐渐这些鸟就会聚成一群，这个群忽高忽低、忽左忽右，直到最后找到食物。这个过程我们转化为一个数学问题。寻找函数　 1(3*x)*()的在[0,4]最大值。该函数的图形如下：

　　当0.9350-0.9450，达到最大值1.3706。为了得到该函数的最大值，我们在[0，4]之间随机的洒一些点，为了演示，我们放置两个点，并且计算这两个点的函数值，同时给这两个点设置在[0，4]之间的一个速度。下面这些点就会按照一定的公式更改自己的位置，到达新位置后，再计算这两个点的值，然后再按照一定的公式更新自己的位置。直到最后在1.3706这个点停止自己的更新。这个过程与粒子群算法作为对照如下：
　　这两个点就是粒子群算法中的粒子。
　　该函数的最大值就是鸟群中的食物　
　　计算两个点函数值就是粒子群算法中的适应值，计算用的函数就是粒子群算法中的适应度函数。
　　更新自己位置的一定公式就是粒子群算法中的位置速度更新公式。
　　下面演示一下这个算法运行一次的大概过程：
　　第一次初始化

第一次更新位置

　　第二次更新位置

　　第21次更新

　　最后的结果（30次迭代）

　　最后所有的点都集中在最大值的地方。
粒子群算法(2)标准的粒子群算法
在上一节的叙述中，唯一没有给大家介绍的就是函数的这些随机的点（粒子）是如何运动的，只是说按照一定的公式更新。这个公式就是粒子群算法中的位置速度更新公式。下面就介绍这个公式是什么。在上一节中我们求取函数1(3*x)*()的在[0,4]最大值。并在[0,4]之间放置了两个随机的点，这些点的坐标假设为x1=1.5； x2=2.5；这里的点是一个标量，但是我们经常遇到的问题可能是更一般的情况－－x为一个矢量的情况，比如二维的情况　2*x1+3*x22的情况。这个时候我们的每个粒子为二维，记粒子P1＝(x1112)2=(x2122)3=(x3132)，(12)。这里n为粒子群群体的规模，也就是这个群中粒子的个数，每个粒子的维数为2。更一般的是粒子的维数为q，这样在这个种群中有n个粒子，每个粒子为q 维。
　　 由n个粒子组成的群体对Q维（就是每个粒子的维数）空间进行搜索。每个粒子表示为：＝（123），每个粒子对应的速度可以表示为(123)，每个粒子在搜索时要考虑两个因素：
　　1。自己搜索到的历史最优值 ，(12)，1,2,3。
　　2。全部粒子搜索到的最优值，(12)，注意这里的只有一个。
　　下面给出粒子群算法的位置速度更新公式：

　　 这里有几个重要的参数需要大家记忆，因为在以后的讲解中将会经常用到：
　　它们是：
是保持原来速度的系数，所以叫做惯性权重。
是粒子跟踪自己历史最优值的权重系数，它表示粒子自身的认识，所以叫“认知”。通常设置为2。
是粒子跟踪群体最优值的权重系数，它表示粒子对整个群体知识的认识，所以叫做“社会知识”，经常叫做“社会”。通常设置为2。
是[0,1]区间内均匀分布的随机数。
是对位置更新的时候，在速度前面加的一个系数，这个系数我们叫做约束因子。通常设置为1。
　　 这样一个标准的粒子群算法就结束了。
　　下面对整个基本的粒子群的过程给一个简单的图形表示：

　　判断终止条件可是设置适应值到达一定的数值或者循环一定的次数。
　　 注意：这里的粒子是同时跟踪自己的历史最优值与全局（群体）最优值来改变自己的位置预速度的，所以又叫做全局版本的标准粒子群优化算法。

粒子群算法(3)标准的粒子群算法(局部版本)
　在全局版的标准粒子群算法中，每个粒子的速度的更新是根据两个因素来变化的，这两个因素是：1. 粒子自己历史最优值。2. 　粒子群体的全局最优值。如果改变粒子速度更新公式，让每个粒子的速度的更新根据以下两个因素更新，A. 粒子自己历史最优值。B. 粒子邻域内粒子的最优值。其余保持跟全局版的标准粒子群算法一样，这个算法就变为局部版的粒子群算法。
　　 一般一个粒子i 的邻域随着迭代次数的增加而逐渐增加，开始第一次迭代，它的邻域为0，随着迭代次数邻域线性变大，最后邻域扩展到整个粒子群，这时就变成全局版本的粒子群算法了。经过实践证明：全局版本的粒子群算法收敛速度快，但是容易陷入局部最优。局部版本的粒子群算法收敛速度慢，但是很难陷入局部最优。现在的粒子群算法大都在收敛速度与摆脱局部最优这两个方面下功夫。其实这两个方面是矛盾的。看如何更好的折中了。
　　 根据取邻域的方式的不同，局部版本的粒子群算法有很多不同的实现方法。
　　第一种方法：按照粒子的编号取粒子的邻域，取法有四种：1，环形取法 2，随机环形取法 3，轮形取法 4，随机轮形取法。
　　 1　 环形2 随机环形
　　 3 轮形 4随机轮形
　　因为后面有以环形取法实现的算法，对环形取法在这里做一点点说明：以粒子1为例，当邻域是0的时候，邻域是它本身，当邻域是1时，邻域为2，8；当邻域是2时，邻域是2，3，7，8；......，以此类推，一直到邻域为4，这个时候，邻域扩展到整个例子群体。据文献介绍（国外的文献），采用轮形拓扑结构，的效果很好。
　　第二种方法：按照粒子的欧式距离取粒子的邻域
　　 　在第一种方法中，按照粒子的编号来得到粒子的邻域，但是这些粒子其实可能在实际位置上并不相邻，于是提出基于空间距离的划分方案，在迭代中计算每一个粒子与群中其他粒子的距离。记录任何2个粒子间的的最大距离为。对每一粒子按照计算一个比值。其中是当前粒子a到b的距离。而选择阈值根据迭代次数而变化。当另一粒子b满足<时，认为b成为当前粒子的邻域。
　　 这种办法经过实验，取得较好的应用效果，但是由于要计算所有粒子之间的距离，计算量大，且需要很大的存储空间，所以，该方法一般不经常使用。

粒子群算法(5)标准粒子群算法的实现
标准粒子群算法的实现思想基本按照粒子群算法(2)标准的粒子群算法的讲述实现。主要分为3个函数。第一个函数为粒子群初始化函数
()其主要作用是初始化粒子群的粒子，并设定粒子的速度、位置在一定的范围内。本函数所采用的数据结构如下所示：
表记录的是粒子的位置、速度与当前的适应度值，我们用W来表示位置，用V来代表速度，用F来代表当前的适应度值。在这里我们假设粒子个数为N，每个粒子的维数为D。
W1,1
W1,2
...
W1
V1,1
V1,2
...
V11
V1
F1
第1个粒子
W2,1
W2,2
...
W2
V2,1
V2,2
...
V21
V2
F2
第2个粒子
...
...
...
...
...

...
...
...
...
.......
1,1
1,2
...
11
1,1
1,2
...
11
1
1
第1个粒子
,1
,2
...

,1
,2
...
1

第N个粒子
表记录每个粒子的历史最优解（粒子历史最好的适应度）以及全部粒子搜索到的全局最优解。用代表全局最优解，W.,1代表每个粒子的历史最优解。粒子群初始化阶段表的前N行与表中的相同，而的值为表中适应度值的最大值对应的行。
,1
,2
...
1

第1个粒子的历史最优解
,1
,2
...
1

第2个粒子的历史最优解
...
...
...
...
...
...
,1
,2
...
1

第1个粒子的历史最优解
,1
,2
...
1

第N个粒子的历史最优解
,1
,2
...
1

全局粒子的历史最优解

根据这样的思想代码如下：
[]()
%功能描述：初始化粒子群，限定粒子群的位置以及速度在指定的范围内
%[]()
%
%输入参数：:种群大小的个数
%输入参数：：一个粒子的维数
%输入参数：:一个粒子在运算中各维的范围；
%　　　　　　　　 格式:
%　　　　　　　　　　 3维粒子的格式:
%　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 [x11
%　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 x22
%　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 x33]
%
%输入参数：：适应度函数
%
%输出：初始化的粒子群
%输出：粒子群当前最优解与全局最优解
%
%用法[]();
%
%异常：首先保证该文件在的搜索路径中，然后查看相关的提示信息。
%
%编制人：
%编制时间：2007.3.26
%参考文献：无
%

%容错控制
4
　　　 ('输入的参数个数错误。')

<2
　　　 ('输出的参数的个数太少，不能保证以后的运行。');

[]();
>1>1
　　　 ('输入的粒子的维数错误，是一个1行1列的数据。');

[]();
2
　　　 ('输入的粒子的维数范围错误。');

%初始化粒子群矩阵

%初始化粒子群矩阵，全部设为[0-1]随机数
('',0);
(,2*1);

%对粒子群中位置,速度的范围进行调节
1
　　　 ()()*((k,2)(k,1))(k,1);
　　　 %调节速度，使速度与位置的范围一致
　　　 ()()*((k,2)(k,1))(k,1);

%对每一个粒子计算其适应度函数的值

1
　　　 (k,2*1)((k,1));

%初始化粒子群最优解矩阵
(1);
%粒子群最优解矩阵全部设为零
[]((:,2*1));
%寻找适应度函数值最大的解在矩阵中的位置(行数)
(1,1);
(1,:)(,1);

下面的函数实现了标准全局版粒子群算法的单步更新位置速度的功能
[]()
%功能描述：全局版本：基本的粒子群算法的单步更新位置,速度的算法
%
%[]()
%
%输入参数：:粒子群矩阵，包含粒子的位置，速度与当前的目标函数值
%输入参数：：包含粒子群个体最优解与全局最优解的矩阵
%输入参数：:一个粒子在运算中各维的范围；
%输入参数：：适应度函数
%输入参数：：迭代的总次数
%输入参数：：当前迭代的次数
%返回值：含意同输入的同名参数
%
%用法：[]()
%
%异常：首先保证该文件在的搜索路径中，然后查看相关的提示信息。
%
%编制人：
%编制时间：2007.3.26
%参考文献：
%参考文献：
%
%修改记录

%2007.3.27
%修改人：
% 添加2*(0,1).*1(,:)中的(0,1)随机数，使性能大为提高
%参照基于的粒子群优化算法程序设计
%
% 总体评价：使用这个版本的调节系数，效果比较好
%

%容错控制
8
　　　 ('输入的参数个数错误。')

2
　　　 ('输出的个数太少，不能保证循环迭代。')

%开始单步更新的操作

%*********************************************
%*****更改下面的代码，可以更改惯性因子的变化*****

%线形递减策略
*(());

固定不变策略
0.7;

%参考文献：陈贵敏，贾建援，韩琪，粒子群优化算法的惯性权值递减策略研究，西安交通大学学报，2006，1
非线形递减，以凹函数递减
()*()^2+()*(2*);

非线形递减，以凹函数递减
*()^(1/(1+10*));
%*****更改上面的代码，可以更改惯性因子的变化*****
%*********************************************

%得到粒子群群体大小以及一个粒子维数的信息
[]();
%得到粒子的维数
(1)/2;
1(1,:)(:,1);

%*********************************************
%*****更改下面的代码，可以更改c12的变化*****
c1=2;
c2=2;

1;
1=4(*(((:,2*1))((1,:))));
2=41;

%*****更改上面的代码，可以更改c12的变化*****
%*********************************************
1
　　 2(1,:)(,1);
　　 .*(1:2*)+2*(0,1).*1(,:)+2*(0,1).*2;
　　 %限制速度的代码
　　 1
　　　　　　 ()>(h,2)
　　　　　　　　　　 ()(h,2);

　　　　　　 ()<(h,2)
　　　　　　　　　　 ()(h,2)+110;
　　　　　　　　　　 %加110防止适应度函数被零除

　　 %更新速度
　　 1:2*);

　　 %*********************************************
　　 %*****更改下面的代码，可以更改约束因子的变化*****

　　 1;

　　 0.729;
　　 %*****更改上面的代码，可以更改约束因子的变化*****
　　 %*********************************************

　　 %限制位置的范围
　　 (,1)*;
　　 1
　　　　　　 ()>(h,2)
　　　　　　　　　　 ()(h,2);

　　　　　　 ()<(h,1)
　　　　　　　　　　 ()(h,1)+110;　　　　　　　　　　

　　 %更新位置
　　 (,1);

　　 %计算每个粒子的新的适应度值
　　 (,2*1)((,1));
　　 (,2*1)>((,1))
　　　　　　 (,1)(,1);

循环结束

%寻找适应度函数值最大的解在矩阵中的位置(行数)，进行全局最优的改变
[]((:,2*1));
((,1))>((1,:))
　　　 (1,:)(,1);

这两个函数给出以后，需要一个函数来把这两个函数组装起来，以此实现一个完整的粒子群算法，这个函数就是
　　代码如下：
[]()
%功能描述：一个循环n次的算法完整过程，返回这次运行的最小与最大的平均适应度,以及在线性能与离线性能
%[]()
%输入参数：:种群大小的个数
%输入参数：：一个粒子的维数
%输入参数：:一个粒子在运算中各维的范围；
%　　　　　　　　 格式:
%　　　　　　　　　　 3维粒子的格式:
%　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 [x11
%　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 x22
%　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 x33]
%
%输入参数:初始化粒子群函数
%输入参数:单步更新速度，位置函数
%输入参数：：适应度函数
%输入参数：：是否每次迭代暂停；＝0，不暂停，否则暂停。缺省不暂停
%输入参数：：是否图形化迭代过程；＝0，不图形化迭代过程，否则，图形化表示。缺省不图形化表示
%输入参数：：迭代的次数；缺省迭代100次
%输入参数：：控制是否绘制在线性能与离线性能的图形表示；0,不显示；
%　　　　　　　　　　　　　　　　 1；显示图形结果。缺省1
%
%返回值：为经过迭代后得到的最优解
%返回值：为在线性能的数据
%返回值：为离线性能的数据
%返回值：为本次完整迭代得到的最小与最大的平均适应度
%
%用法[]();
%
%异常：首先保证该文件在的搜索路径中，然后查看相关的提示信息。
%
%编制人：
%编制时间：2007.3.26
%参考文献：

%修改记录：
%添加，以得到性能评估数据
%修改人：
%修改时间：2007.3.27
%参考文献：.

%容错控制
<4
　　　 ('输入的参数个数错误。')

[]();
>1>1
　　　 ('输入的粒子的维数错误，是一个1行1列的数据。');

[]();
2
　　　 ('输入的粒子的维数范围错误。');

%设置缺省值
<7
　　　 1;
　　　 100;
　　　 0;
　　　 0;

<8
　　　 1;
　　　 0;
　　　 100;

<9
　　　 100;
　　　 1;

<10
　　　 1;

%控制是否显示2维以下粒子维数的寻找最优的过程
0
　　　 ();

%初始化种群　　　　　　
[]()

%在测试函数图形上绘制初始化群的位置
0
　　　 1
　　　 1
　　　　　　　 ([(,1)(,1)],[(,3),0],'r*-','',8);
　　　　　　　 ((,1)(,3)2());

　　　 2
　　　　　　　 1
　　　　　　　　　　　 3((,1)(,2)(,5),'r.','',8);

%暂停让抓图
0
　　　 ('开始迭代，按任意键：')

%开始更新算法的调用
1
　　　 %显示迭代的次数：
　　　 ('')
　　　 ('第 此迭代');
　　　 ();
　　　 ('')

　　　 %调用一步迭代的算法
　　　 [](,0.95,0.4)

　　　 %在目标函数的图形上绘制2维以下的粒子的新位置
　　　 0
　　　　　　　 1
　　　　　　　　　　　 1
　　　　　　　　　　　　　　　 ([(,1)(,1)],[(,3),0],'r*-','',8);
　　　　　　　　　　　　　　　 ((,1)(,3)2());

　　　　　　　 2
　　　　　　　　　　　 1
　　　　　　　　　　　　　　　 3((,1)(,2)(,5),'r.','',8);

　　　 (1,1);
　　　 ();　　　
　　　 0
　　　　　　　 (1,1);
　　　　　　　 ();　　　
　　　　　　　 ('步迭代的最优目标函数值');
　　　　　　　 ();
　　　　　　　 ('下次迭代，按任意键继续');

　　　 %记录每一步的平均适应度
　　　 (1)((:,2*1));

循环结束标志

%记录最小与最大的平均适应度
[()()];
%计算离线与在线性能
1
　　　 (1)((1,1));
　　　 (1)((1,1));

1
　　　 (1)((1,1));

%绘制离线性能与在线性能曲线
1

　　　 ('离线性能曲线图')
　　　 ('迭代次数');
　　　 ('离线性能');

　　　 ();

　　　 ('在线性能曲线图')
　　　 ('迭代次数');
　　　 ('在线性能');

　　　 ();

%记录本次迭代得到的最优结果
(1,1);
();
[];
这里给出一个使用的例子代码，并分别解释各参数的含义：
%打开计时器
;
%
[-50 50
　　　 -50 50
　　　 -50 50
　　　 -50 50
　　　 -50 50
　　　 -50 50
　　　 -50 50
　　　 -50 50
　　　 -50 50
　　　 -50 50];
[](20,10,0,0,4000,0);

在上面的代码中函数中的20代表粒子群的规模为20个，10代表每个粒子的维数为10，是粒子的每一维的范围，同时也是速度的范围，是初始化函数的句柄，是单步更新的函数句柄，是适应度评价函数的句柄，4000代表真个算法循环4000次终止，其他参数参见说明文档。
粒子群算法(6)几个适应度评价函数
　下面给出几个适应度评价函数，并给出图形表示
　　 头几天机子种了病毒，重新安装了系统，不小心把程序全部格式化了，痛哭！！！没办法，好多程序不见了，现在把这几个典型的函数重新编写了，把他们给出来，就算粒子群算法的一个结束吧！痛恨病毒！！！！

　　第一个函数：函数，图形如下所示：

适应度函数如下：（为了求最大值，我去了所有函数值的相反数）　
　(x)
函数
%输入x,给出相应的y值,在(0,0,…,0)处有全局极小点0.
%编制人：
%编制日期：
[](x);
　>1
　　　　('输入的参数错误');

y1=1/4000*(x.^2);
y2=1;
　1
　　　　y22*(x(h)(h));

12+1;
;
绘制函数图像的代码如下：
　()
%绘制函数图形
[-8:0.1:8];
;
[]();
[](X);
　1
　　　　　1
　　　　　　　　z()([X()()]);

();

绘制函数图像的代码如下：

　()
%绘制函数图形
[-5:0.05:5];
;
[]();
[](X);
　1
　　　　　1
　　　　　　　　z()([X()()]);

();