粒子群算法c语言

‘壹’ 优化算法笔记（五）粒子群算法（3）

(已合并本篇内容至粒子群算法（1）)

上一节中，我们看到小鸟们聚集到一个较小的范围内后，不会再继续集中。这是怎么回事呢？
猜测：
1.与最大速度限制有关，权重w只是方便动态修改maxV。
2.与C1和C2有关，这两个权重限制了鸟儿的搜索行为。
还是上一节的实验， 。现在我们将maxV的值有5修改为50，即maxV=50，其他参数不变。参数如下

此时得到的最优位值的适应度函数值为0.25571,可以看出与maxV=5相比，结果差了很多而且小鸟们聚集的范围更大了。
现在我们设置maxV=1，再次重复上面的实验，实验结果如下：

这次最终的适应度函数值为，比之前的结果都要好0.00273。从图中我们可以看出，小鸟们在向一个点集中，但是他们飞行的速度比之前慢多了，如果问题更复杂，可能无法等到它们聚集到一个点，迭代就结束了。
为什么maxV会影响鸟群的搜索结果呢？
我们依然以maxV=50为例，不过这次为了看的更加清晰，我们的鸟群只有2只鸟，同时将帧数放慢5倍以便观察。

思路一：限制鸟的最大飞行速率,由于惯性系数W的存在，使得控制最大速率和控制惯性系数的效果是等价的，取其一即可。
方案1：使惯性系数随着迭代次数增加而降低，这里使用的是线性下降的方式，即在第1次迭代，惯性系数W=1,最后一次迭代时，惯性系数W=0，当然，也可以根据自己的意愿取其他值。
实验参数如下：

小鸟们的飞行过程如上图，可以看到效果很好，最后甚至都聚集到了一个点。再看看最终的适应度函数值8.61666413451519E-17，这已经是一个相当精确的值了，说明这是一个可行的方案，但是由于其最后种群过于集中，有陷入局部最优的风险。
方案2：给每只鸟一个随机的惯性系数，那么鸟的飞行轨迹也将不再像之前会出现周期性。每只鸟的惯性系数W为（0,2）中的随机数（保持W的期望为1）。
实验参数如下：

可以看到小鸟们并没有像上一个实验一样聚集于一个点，而是仍在一个较大的范围内进行搜索。其最终的适应度函数为0.01176，比最初的0.25571稍有提升，但并不显着。什么原因造成了这种情况呢？我想可能是由于惯性系数成了期望为1的随机数，那么小鸟的飞行轨迹的期望可能仍然是绕着一个四边形循环，只不过这个四边形相比之前的平行四边形更加复杂，所以其结果也稍有提升，当然对于概率算法，得到这样的结果可能仅仅是因为运气不好
我们看到惯性系数W值减小，小鸟们聚拢到一处的速度明显提升，那么，如果我们去掉惯性系数这个参数会怎么样呢。
方案3：取出惯性系数，即取W=0，小鸟们只向着那两个最优位置飞行。

可以看见鸟群们迅速聚集到了一个点，再看看得到的结果，最终的适应度函数值为2.9086886073362966E-30，明显优于之前的所有操作。
那么问题来了，为什么粒子群算法需要一个惯性速度，它的作用是什么呢？其实很明显，当鸟群迅速集中到了一个点之后它们就丧失了全局的搜索能力，所有的鸟会迅速向着全局最优点飞去，如果当前的全局最优解是一个局部最优点，那么鸟群将会陷入局部最优。所以，惯性系数和惯性速度的作用是给鸟群提供跳出局部最优的可能性，获得这个跳出局部最优能力的代价是它们的收敛速度减慢，且局部的搜索能力较弱（与当前的惯性速度有关）。
为了平衡局部搜索能力和跳出局部最优能力，我们可以人为的干预一下惯性系数W的大小，结合方案1和方案2，我们可以使每只鸟的惯性系数以一个随机周期，周期性下降，若小于0，则重置为初始值。

这样结合了方案1和方案2的惯性系数，也能得到不错的效果，大家可以自己一试。

思路二：改变小鸟们向群体最优飞行和向历史最优飞行的权重。
方案4：让小鸟向全局最优飞行的系数C2线性递减。

小鸟们的飞行过程与之前好像没什么变化，我甚至怀疑我做了假实验。看看最终结果，0.7267249621552874，这是到目前为止的最差结果。看来这不是一个好方案，让全局学习因子C2递减，势必会降低算法的收敛效率，而惯性系数还是那么大，小鸟们依然会围绕历史最优位置打转，毕竟这两个最优位置是有一定关联的。所以让C1线性递减的实验也不必做了，其效果应该与方案4相差不大。
看来只要是惯性系数不变怎么修改C1和C2都不会有太过明显的效果。为什么实验都是参数递减，却没有参数递增的实验呢？
1.惯性系数W必须递减，因为它会影响鸟群的搜索范围。
2.如果C1和C2递增，那么小鸟的惯性速度V势必会跟着递增，这与W递增会产生相同的效果。

上面我们通过一些实验及理论分析了粒子群算法的特点及其参数的作用。粒子群作为优化算法中模型最简单的算法，通过修改这几个简单的参数也能够改变算法的优化性能可以说是一个非常优秀的算法。
上述实验中，我们仅分析了单个参数对算法的影响，实际使用时（创新、发明、写论文时）也会同时动态改变多个参数，甚至是参数之间产生关联。
实验中，为了展现实验效果，maxV取值较大，一般取值为搜索空间范围的10%-20%，按上面（-100,100）的范围maxV应该取值为20-40，在此基础上，方案1、方案2效果应该会更好。
粒子群算法是一种概率算法，所以并不能使用一次实验结果来判断算法的性能，我们需要进行多次实验，然后看看这些实验的效果最终来判断，结果必须使用多次实验的统计数据来说明，一般我们都会重复实验30-50次，为了发论文去做实验的小伙伴们不要偷懒哦。
粒子群算法的学习目前告一段落，如果有什么新的发现，后面继续更新哦！
以下指标纯属个人yy,仅供参考

目录
上一篇 优化算法笔记（四）粒子群算法（2）
下一篇 优化算法笔记（六）遗传算法

‘贰’ C语言，用粒子群算法求函数最大值，拜托了！！

for i=1:sizepop % 随机产生一个种群 pop(i,:)=2*rands(1,2); % 初始化粒子 V(i,:)=0.5*rands(1,2); % 初始化速度 % 计算粒子适应度值 fitness(i)=fun(pop(i,:)); end [bestfitness bestindex]=min(fitness); zbest=pop(bestindex,:); % 群体极...

‘叁’ 粒子群算法简单介绍

粒子群算法（也称粒子群优化算法（particle swarm optimization, PSO）），模拟鸟群随机搜索食物的行为。粒子群算法中，每个优化问题的潜在解都是搜穗穗索空间中的一只鸟，叫做“粒子”。所有的粒子都有一个由被优化的函数决定的适应值（fitness value），每个粒子还有一个速度决定它们“飞行”的方向和距离。

粒子群算法初始化为一群随机的粒子（随机解），然后根据迭代找到最优解。每一次迭代中，粒子通过跟踪两个极敏扰值来更新自己：第1个是粒子本身所找到的最优解，这个称为个体极值；第2个是整个种群目前找到的最优解，这个称为全局极值。也可以不用整个种群，而是用其中的一部分作为粒子的邻居猜拿卜，称为局部极值。

假设在一个D维搜索空间中，有N个粒子组成一个群落，其中第i个粒子表示为一个D维的向量：

第i个粒子的速度表示为：

还要保存每个个体的已经找到的最优解 ，和一个整个群落找到的最优解 。

第i个粒子根据下面的公式更新自己的速度和位置：

其中， 是个体已知最优解， 是种群已知最优解， 为惯性权重， , 为学习因子（或加速常数 acceleration constant)， , 是[0,1]范围内的随机数。

式(1)由三部分组成：

‘肆’ 粒子群算法

粒子群算法（Particle Swarm Optimization），又称鸟群觅食算法，是由数学家J. Kennedy和R. C. Eberhart等开发出的一种新的进化算法。它是从随机解开始触发，通过迭代寻找出其中的最优解。本算法主要是通过适应度来评价解的分数，比传统的遗传算法更加的简单，它没有传统遗传算法中的“交叉”和“变异”等操作，它主要是追随当前搜索到的最优值来寻找到全局最优值。这种算法实现容易，精度高，收敛快等特点被广泛运用在各个问题中。

粒子群算法是模拟鸟群觅食的所建立起来的一种智能算法，一开始所有的鸟都不知道食物在哪里，它们通过找到离食物最近的鸟的周围，再去寻找食物，这样不断的追踪，大量的鸟都堆积在食物附近这样找到食物的几率就大大增加了。粒子群就是这样一种模拟鸟群觅食的过程，粒子群把鸟看成一个个粒子，它们拥有两个属性——位置和速度，然后根据自己的这两个属性共享到整个集群中，其他粒子改变飞行方向去找到最近的区域，然后整个集群都聚集在最优解附近，最后最终找到最优解。

算法中我们需要的数据结构，我们需要一个值来存储每个粒子搜索到的最优解，用一个值来存储整个群体在一次迭代中搜索到的最优解，这样我们的粒子速度和位置的更新公式如下：

其中pbest是每个粒子搜索到的最优解，gbest是整个群体在一次迭代中搜索到的最优解，v[i]是代表第i个粒子的速度，w代表惯性系数是一个超参数，rang()表示的是在0到1的随机数。Present[i]代表第i个粒子当前的位置。我们通过上面的公式不停的迭代粒子群的状态，最终得到全局最优解

‘伍’ 关于粒子群算法的问题

粒子群的版本甚多，常用的是加有惯性权重w的
v[] = w * v[] + c1 * rand() * (pbest[] - present[]) + c2 * rand() * (gbest[] - present[])
一般选择惯性权重在迭代过程中线性下降，目的是在迭代的初期，以比较大的权重分配给粒子的原速度，而防止粒子过早的倾向于其本身的局部最优与全局最优，此时的全局搜索能力是可以的。但粒子群是基于牛顿力学的，随着w的减小，速度v的作用会在更新中弱化，对应的是，pbest和gbest的作用得到了加强，这也就意味着，粒子会更加趋向于pbest和gbest的方向移动。这个时候粒子就特别容易陷入局部最优了。

其实陷入局部最优不只是粒子群的问题，进化类的算法都存在这个问题，只不过有些算法随机性强一些，收敛速度慢一些，所以更加容易跳出局部最优（但不是绝对避免）

‘陆’ 粒子群算法的算法介绍

如前所述，PSO模拟鸟群的捕食行为。设想这样一个场景：一群鸟在随机搜索食物。在这个区域里只有一块食物。所有的鸟都不知道食物在那里。但是他们知道当前的位置离食物还有多远。那么找到食物的最优策略是什么呢。最简单有效的就是搜寻目前离食物最近的鸟的周围区域。
PSO从这种模型中得到启示并用于解决优化问题。PSO中，每个优化问题的解都是搜索空间中的一只鸟。我们称之为“粒子”。所有的粒子都有一个由被优化的函数决定的适应值(fitness value)，每个粒子还有一个速度决定他们飞翔的方向和距离。然后粒子们就追随当前的最优粒子在解空间中搜索。
PSO 初始化为一群随机粒子(随机解)。然后通过迭代找到最优解。在每一次迭代中，粒子通过跟踪两个极值来更新自己。第一个就是粒子本身所找到的最优解，这个解叫做个体极值pBest。另一个极值是整个种群目前找到的最优解，这个极值是全局极值gBest。另外也可以不用整个种群而只是用其中一部分作为粒子的邻居，那么在所有邻居中的极值就是局部极值。 在找到这两个最优值时，粒子根据如下的公式来更新自己的速度和新的位置：
v[] = w * v[] + c1 * rand() * (pbest[] - present[]) + c2 * rand() * (gbest[] - present[]) (a)
present[] = present[] + v[] (b)
v[] 是粒子的速度, w是惯性权重,present[] 是当前粒子的位置. pbest[] and gbest[] 如前定义 rand () 是介于（0， 1）之间的随机数. c1, c2 是学习因子. 通常 c1 = c2 = 2.
程序的伪代码如下
For each particle
____Initialize particle
END
Do
____For each particle
________Calculate fitness value
________If the fitness value is better than the best fitness value (pBest) in history
____________set current value as the new pBest
____End
____Choose the particle with the best fitness value of all the particles as the gBest
____For each particle
________Calculate particle velocity according equation (a)
________Update particle position according equation (b)
____End
While maximum iterations or minimum error criteria is not attained
在每一维粒子的速度都会被限制在一个最大速度Vmax，如果某一维更新后的速度超过用户设定的Vmax，那么这一维的速度就被限定为Vmax

‘柒’ 求粒子群算法C#版

static void Main(string[] args)
{
PSO p = new PSO();
p.mains();
System.Console.ReadLine();
}
class PSO
{
const int S = 20; /*试验次数 */
const int G = 2000; /*混合迭代次数*/
const int P = 40; /*个体总数*/
//#define c1 2.0 /*学习因子1*/
//#define c2 2.0 /*学习因子2*/
const int V = 30; /*个体维数*/
const double MAX = 5.12;
const double MIN = -5.12;
double D = MAX; /*蛙跳的最大值*/
int i1, i2, i3, i4;
int try_number = 0;
int try_max = 5;
double R;//0-1之间的随机数,精度为1/10000
double Wmax = 0.9;
double Wmin = 0.4;
double PI = 3.14159265;
double W = 0.9;
double c1 = 2.0;/*学习因子1*/
double c2 = 2.0;/*学习因子2*/
double Vpso = 0.0;
double Tolerance = 0.0000001;//收敛精度
double c3 = 0.03;//扰动幅度
double e = 2.718281828459;//自然对数底数
int sm = 3;
int bz = 0;//扰动因子标志
public class Individal
{
public double[] d = new double[V];
public double fitness;
}
public class psom
{
public double[] sd = new double[V];
} public psom[] pso = new psom[P]; /*记录每个离子各位的更新速度*/
public Individal px; /*全体中最好位置*/
public Individal[] indivial = new Individal[P]; /*全部个体*/
public Individal[] indiviala = new Individal[P]; /*全部个体——备份*/
public Individal tem; public PSO()
{
for (int i = 0; i < P; i++)
{
pso[i] = new psom();
indivial[i] = new Individal();
indiviala[i] = new Individal();
}
}
/*选择测试函数为Sphere*/
/*选择测试函数为Sphere*/
public double fitness(double[] a)
{
int i;
double sum = 0.0;
double sum1 = 0.0;
double s1 = 0.0, h1 = 0.0;
double[] x1 = new double[V + 1];
for (i = 0; i < V; i++) x1[i] = a[i];
for (i = 0; i < V; i++)
for (i = 0; i < V; i++)
sum = sum + (x1[i] * x1[i] - 10 * Math.Cos(2 * PI * x1[i]) + 10); return sum;
} /*对每一个个体初始化*/
public void init()
{
Random ran = new Random();
//R = ran.NextDouble(); int i, j, pmin = 0;
//srand((unsigned)time(NULL));
for (i = 0; i < P; i++)
{
for (j = 0; j < V; j++)
{
R = ran.NextDouble();
indivial[i].d[j] = R * (MAX - MIN) + MIN;
}
indivial[i].fitness = fitness(indivial[i].d);//计算初始适应值
indiviala[i] = indivial[i];//将个体复制给另一个序列
}
for (i = 0; i < P - 1; i++)
{
if (indivial[pmin].fitness > indivial[i + 1].fitness)
pmin = i + 1;//适应值小者最优
}
px = indivial[pmin];//最优者为小
}
/*按照适应度降序对全部个体进行排序和族群划分*/ /*群组内更新*/
public void update()
{
int i, j, k, l, n;
double a;
double b;
W = Wmax - (double)(i2) * (Wmax - Wmin) / (double)(G);
for (i = 0; i < P; i++)
{ for (j = 0; j < V; j++)//更新粒子速度、位置
{//更新速度
pso[i].sd[j] = W * pso[i].sd[j] + c1 * R * (indiviala[i].d[j] - indivial[i].d[j]) + c2 * R * (px.d[j] - indivial[i].d[j]);
if (pso[i].sd[j] > D) pso[i].sd[j] = D; //D 最大速度
if (pso[i].sd[j] < -D) pso[i].sd[j] = -D;
indivial[i].d[j] = indivial[i].d[j] + pso[i].sd[j];//更新位置
}
a = fitness(indivial[i].d);//计算本次迭代的粒子适应值
// printf("old是%.16f",indivial[i].fitness);
indivial[i].fitness = a;
// printf("new是%.16f",indivial[i].fitness);
// getchar();
if (a < indiviala[i].fitness)
{
indiviala[i] = indivial[i];
if (indiviala[i].fitness < px.fitness)
px = indiviala[i];
}//比较粒子与前一次迭代的适应值 寻求最优者
}
}
public void report()
{
int i;
System.Console.WriteLine(px.fitness);

for (i = 0; i < P; i++)
{
// printf("%.16f\n",indivial[i].fitness);
// printf("%.16f\n",indiviala[i].fitness);
}
}
public void mains()
{
// int i1,i2;
//clock_t start, end;
double ave;
//FILE* f = fopen("result(SFLA).txt", "w");
//for(i4=0;i4<G+1;i4+=50)
//{
ave = 0.0;
//start = clock();
for (int i1 = 0; i1 < S; i1++)
{ init();
for (int i2 = 0; i2 < G; i2++)
{
update();
// report();
// getchar();
}
//
//
report();
ave = ave + px.fitness;
}
//end = clock();
ave = ave / S;
System.Console.WriteLine("平均极值为:");
System.Console.WriteLine(ave);
//System.Console.WriteLine("Interval=%.2fseconds\n", (double)(end - start) / ((double)CLOCKS_PER_SEC));
//printf("%d代为平均极值为%.16f\n",i4,ave);
//fprintf(f,"%d代为平均极值为%.16f\n",i4,ave);
// getchar();
//getchar();
}
}

‘捌’ 粒子群算法（一）：粒子群算法概述

  本系列文章主要针对粒子群算法进行介绍和运用，并给出粒子群算法的经典案例，从而进一步加深对粒子群算法的了解与运用（预计在一周内完成本系列文章）。主要包括四个部分：

  粒子群算法也称粒子群优化算法（Particle Swarm Optimization, PSO），属于群体智能优化算法，是近年来发展起来的一种新的进化算法（Evolutionary Algorithm, EA）。 群体智能优化算法主要模拟了昆虫、兽群、鸟群和鱼群的群集行为，这些群体按照一种合作的方式寻找食物，群体中的每个成员通过学习它自身的经验和其他成员的经验来不断地改变搜索的方向。 群体智能优化算法的突出特点就是利用了种群的群体智慧进行协同搜索，从而在解空间内找到最优解。
  PSO 算法和模拟退火算法相比，也是 从随机解出发，通过迭代寻找最优解 。它是通过适应度来评价解的品质，但比遗传算法规则更为简单，没有遗传算法的“交叉”和“变异”，它通过追随当前搜索到的最大适应度来寻找全局最优。这种算法以其 容易实现、精度高、收敛快 等优点引起了学术界的重视，并在解决实际问题中展示了其优越性。

  在粒子群算法中，每个优化问题的解被看作搜索空间的一只鸟，即“粒子”。算法开始时首先生成初始解，即在可行解空间中随机初始化 粒子组成的种群 ，其中每个粒子所处的位置 ，都表示问题的一个解，并依据目标函数计算搜索新解。在每次迭代时，粒子将跟踪两个“极值”来更新自己， 一个是粒子本身搜索到的最好解 ，另一个是整个种群目前搜索到的最优解 。 此外每个粒子都有一个速度 ，当两个最优解都找到后，每个粒子根据如下迭代式更新：

  其中参数 称为是 PSO 的 惯性权重(inertia weight) ，它的取值介于[0,1]区间；参数 和 称为是 学习因子(learn factor) ；而 和 为介于[0,1]之间的随机概率值。
  实践证明没有绝对最优的参数，针对不同的问题选取合适的参数才能获得更好的收敛速度和鲁棒性，一般情况下 , 取 1.4961 ，而 采用 自适应的取值方法 ，即一开始令 ， 使得 PSO 全局优化能力较强 ；随着迭代的深入，递减至 ， 从而使得PSO具有较强的局部优化能力 。

  参数 之所以被称之为惯性权重，是因为 实际 反映了粒子过去的运动状态对当前行为的影响，就像是我们物理中提到的惯性。 如果 ，从前的运动状态很少能影响当前的行为，粒子的速度会很快的改变；相反， 较大，虽然会有很大的搜索空间，但是粒子很难改变其运动方向，很难向较优位置收敛，由于算法速度的因素，在实际运用中很少这样设置。也就是说， 较高的 设置促进全局搜索，较低的 设置促进快速的局部搜索。