粒子群算法

算法介绍

简介

如前所述，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.

C#源代码如下：

using System;using System.Linq;namespace MSAlgorithm{    /// <summary>    /// 粒子群算法,也称粒子群优化算法（Particle Swarm Optimization）    /// </summary>    public class PSO    {        /// <summary>        /// 粒子数        /// </summary>        const int NUM = 40;        /// <summary>        /// 维数        /// </summary>        const int DIM = 30;        /// <summary>        /// 参数        /// </summary>        const double C1 = 1.8;        /// <summary>        /// 参数        /// </summary>        const double C2 = 1.8;        /// <summary>        /// 位置下限        /// </summary>        static double xMin = -100.0;        /// <summary>        /// 位置上限        /// </summary>        static double xMax = 100.0;        /// <summary>        /// 全局最优位置        /// </summary>        static double[] gbestX = new double[DIM];        /// <summary>        /// 全局最优适应度        /// </summary>        static double gbestF;        /// <summary>        /// 用于生成随机数        /// </summary>        static Random rand = new Random();        /// <summary>        /// 定义一个粒子类        /// </summary>        class Particle        {            /// <summary>            /// 当前位置矢量            /// </summary>            public double[] x = new double[DIM];            /// <summary>            /// 历史最优位置            /// </summary>            public double[] bestX = new double[DIM];            /// <summary>            /// 当前适应度            /// </summary>            public double f;            /// <summary>            /// 历史最优适应度            /// </summary>            public double bestF;        }        /// <summary>        /// 内部函数,超球函数        /// </summary>        /// <param name="x">当前位置矢量</param>        /// <returns>位置</returns>        static double f1(double[] x)        {            return x.Sum(a => a * a);        }        /// <summary>        /// 获取全局最优适应度        /// </summary>        /// <returns></returns>        public double GetGbestF()        {            //定义粒子群            Particle[] swarn = new Particle[NUM];            //初始化全局最优            for (int i = 0; i < DIM; i++)                gbestX[i] = rand.NextDouble() * (xMax - xMin) + xMin;            gbestF = double.MaxValue;            //初始化粒子群            for (int j = 0; j < NUM; j++)            {                Particle p1 = new Particle();                for (int d = 0; d < DIM; d++)                    p1.x[d] = rand.NextDouble() * (xMax - xMin) + xMin;                p1.f = f1(p1.x);                p1.bestF = double.MaxValue;                swarn[j] = p1;            }            for (int t = 0; t < 5000; t++)            {                for (int n = 0; n < NUM; n++)                {                    Particle p1 = swarn[n];                    //进化方程                    for (int d = 0; d < DIM; d++)                        p1.x[d] += C1 * rand.NextDouble() * (p1.bestX[d] - p1.x[d]) + C2 * rand.NextDouble() * (gbestX[d] - p1.x[d]);                    p1.f = f1(p1.x);                    //改变历史最优                    if (p1.f < p1.bestF)                    {                        p1.x.CopyTo(p1.bestX, 0);                        p1.bestF = p1.f;                    }                    //改变全局最优                    if (p1.f < gbestF)                    {                        p1.x.CopyTo(gbestX, 0);                        //把当前全局最优的粒子随机放到另一位置                        for (int d = 0; d < DIM; d++)                            p1.x[d] = rand.NextDouble() * (xMax - xMin) + xMin;                        gbestF = p1.f;                    }                }            }            return gbestF;        }    }}

运行的结果图示：


参考文献：
百度百科