粒子群算法

算法介绍

简介

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

运行的结果图示：


参考文献：
百度百科