用粒子群优化算法求解旅行商问题

  • 演示程序下载 - 116.2 KB

前言

粒子群优化算法采用一种人工智能的形式来解决问题。这种算法对于求解那些使用了多个连续变化的值的函数来说,尤为有效。这篇文章将会介绍如何修改粒子群算法,以使用离散固定值来解决诸如旅行商(TSP,Travelling Salesman Problem)这样的问题。

背景知识

关于粒子群优化算法(PSO,Particle Swarm Optimizers),我在以前的文章中已经进行过讨论与论证。正如我在那篇文章中所说,这一算法的基本思想,是在问题所有可能的解决方案的范围内移动(飞行)解决问题的实体(粒子)的群(群)。我们将这一范围称作问题空间。粒子在问题空间内的运动具有随机性,而这随机性则主要取决于以下三个因素:

  1. 粒子的当前所处位置。
  2. 粒子所找到的最好的位置,称为“个体最优”(Personal Best)或用 pBest 表示。
  3. 在群中发现的最好的位置,则称为“全局最优”(Global Best)或用 gBest 表示。

该算法的现代变体一般使用的是局部最佳位置而非全局最佳位置,因为这往往能够更好地探索问题空间,防止过快地收敛到某些区域的最小值。在这些变体中,群体被分成一些被称为“线人(Informers)”的粒子群。信息在小组成员内部互相交换,从而可以确定该小组的局部最优位置。当出现一定迭代次数后,全局最优值没有发生改变的情况,则要将粒子重组以获得新的组合。

最初的粒子群优化算法公式

处理连续变量的公式为:

**·** **Vid = vid \* W + C1 \* rand(pid - xid)+ C2 \* Rand(pgd - xid)**

其中,vid是当前速度,Vid是新速度。在我们所描述的情况中,速度即是位置改变的量。

W,C1,C2 是常数。常数的估计值为 C1 = C2 = 1.4,W = 0.7。

Randrand 则是两个随机生成的数值,其值域都为 [0, 1)。

xid 是当前位置,pid 是个体最优位置,pgd 是全局最优位置。我们通过对位置添加新的速度来更新它。xid = xid + Vid,此公式适用于位置的每个维度。

旅行商问题简介

如何找到一个总距离最短的行走路径,这一路径使得旅行商访问了每一个城市,且每个城市仅访问了一次,最后还要回到他最开始的位置。这并非一个完全学术化的问题,在接线图和印刷电路板的设计中我们也会碰到类似的情况。这是一个有着许多标准的[城市列表](http://elib.zib.de/pub/mp-testdata/tsp/tsplib/tsp/)例子的,记录完备的问题。当下已经有许多关于如何使用 PSO 来解决这个问题的论文。而我们这里所使用的方法是基于一篇名为《[结合遗传算法与粒子群算法解决旅行商问题](https://www.cogentoa.com/article/10.1080/23311835.2015.1048581.pdf)》的文章,作者是 Keivan Borna 和 Razieh Khezri。

使用 PSO 更新旅行商的路线

正如我们所看到的,粒子的新位置会受到三个因素的不同程度的影响。这三个因素分别是粒子当前位置,前一个最优位置,以及组内的最优位置。旅行商的路线可以根据这些因素,分成三个部分进行更新,其中每个部分的大小由该部分的相对强度决定。这些部分可以连接在一起形成一个新路线。但这种方法存在一个问题,因为每个城市只能被列入一次,然而某部分路径可能包含了已经在之前的路线中列出的城市。所以我们需要有一个确保每个城市都被加入到这个路线中,并且在这个过程中没有任一城市重复的机制。

在上图中,我们从当前路径选择的部分是 6, 3, 5。则这些城市会被添加到新的路径。个人最优路径中选择了城市 1, 3, 2,由于城市 3 已被添加,所以我们只选择城市 1 和 2。局部最优路径选定了城市 7, 3,当然由于城市 3 已经被添加,所以只有城市 7 被选中。最后,未被选中的两个城市,即城市 0 和城市 4 ,就按照它们出现在当前路径中的顺序添加到新路径当中。

通过使用 [BitArrays](https://msdn.microsoft.com/en-us/library/system.collections.bitarray(v=vs.110%29),我们所选择的要添加的城市很方便就能加入数组中。一个`BitArray`用作可用性掩码,它所有的位最初设置为 true。而另一个`BitArray`则用于对需要添加的片段的选择掩码。为了说明这一点,请看下面这个添加当前分段后的情况。

示例程序

随机数生成

该应用程序需要产生大量的随机数,所以为此寻找一个最佳的随机数生成器(RNG,Random Number Generator)是值得的。经过大量的研究,我发现`System.Random` 和其它一些好的生成器一样强大,甚至还比大部分方法好。如果您有兴趣研究 RNG 的质量,那么在这里有一个用 C# 编写的 15 个测试的 [Diehard ](https://en.wikipedia.org/wiki/Diehard_tests)序列的链接(译者注:该链接已失效,故此处无超链接)。由于某种原因,我无法运行测试 2,也许是因为我的资源比所需求的 8000 万比特少了点吧。

城际查找表

为找到两个城市之间的距离,程序使用了一张二维矩阵形式的查找表。例如,为了得到城市 A 和城市 B 之间的距离,查找城市 A 对应的行和城市 B 对应的列,在行和列所确定的位置处就能获取距离信息。该表是以 [Indexer ](https://docs.microsoft.com/en-us/dotnet/csharp/programming-guide/indexers/)的形式实现的,因此它实际上就成了一个只读的二维数组。让其只读是因为考虑到表是由多个对象共享的,所以最好令它是不可变的。

public class LookUpTable<T> : ILookUpTable<T>
    {
        public LookUpTable(T[,] arr)
        {
            this.arr = arr;
        }

        private readonly T[,] arr;

        // The indexer allows the use of [,] operator.
        public T this[int r, int c]
        {
            get
            {
                return arr[r, c];
            }
        }

        public int RowCount
        {
            get
            {
                return arr.GetLength(0);
            }
        }
        public int ColCount
        {
            get
            {
                return arr.GetLength(1);
            }
        }
    }

实现

示例应用程序使用`TspParticle`对象数组来实现一个群,并采用`SwarmOptimizer`来优化它。程序将从文件 app.config 中读入优化算法相关的属性,例如群大小(Swarm Size)和循环次数(Epoch)。

public int Optimize(ITspParticle[] swarm, PsoAttributes psoAttributes)
       {
           this.BestGlobalItinery = swarm[0].CurrentRoute.Clone();
           int[] particleIndex = this.InitArray(psoAttributes.SwarmSize);
           int epoch = 0;
           int staticEpochs = 0;
           while (epoch < psoAttributes.MaxEpochs)
           {
               bool isDistanceImproved = false;
               foreach (ITspParticle particle in swarm)
               {
                   double distance = particle.Optimize(psoAttributes);
                   if (distance < this.BestGlobalItinery.TourDistance)
                   {
                       particle.CurrentRoute.CopyTo(this.BestGlobalItinery);
                       isDistanceImproved = true;
                       this.consoleWriter.WriteDistance(distance);
                   }
               }
               if (!isDistanceImproved)
               {
                   staticEpochs++;
                   if (staticEpochs == psoAttributes.MaxStaticEpochs)
                   {
                       this.UpdateInformers(swarm, particleIndex, psoAttributes.MaxInformers);
                       staticEpochs = 0;
                   }
               }
               epoch++;
           }
           return (int)this.BestGlobalItinery.TourDistance;
       }

每个粒子都包含对它的`CurrentRoute`,`PersonalBestRoute`与`LocalBestRoute`的引用,它们以整形数组的形式包含要访问的城市的顺序,此处最后一个列出的城市会链接到开始的第一个城市。路径使用一个`ParticleOptimizer`来进行更新。

public int[] GetOptimizedDestinationIndex(
 IRoute currRoute,
 IRoute pBRoute,
 IRoute lBRoute,
 PsoAttributes psoAttribs)
  {
    //update all the velocities using the appropriate PSO constants
    double currV = routeManager.UpdateVelocity(currRoute, psoAttribs.W, 1);
    double pBV = routeManager.UpdateVelocity(pBRoute, psoAttribs.C1, randomFactory.NextRandomDouble());
    double lBV = routeManager.UpdateVelocity(lBRoute, psoAttribs.C2, randomFactory.NextRandomDouble());
    double totalVelocity = currV + pBV + lBV;

    //update the Segment size for each Route
    currRoute.SegmentSize = routeManager.GetSectionSize(currRoute, currV, totalVelocity);
    pBRoute.SegmentSize = routeManager.GetSectionSize(pBRoute, pBV, totalVelocity);
    lBRoute.SegmentSize = routeManager.GetSectionSize(lBRoute, lBV, totalVelocity);
    return routeManager.AddSections(new[] { lBRoute, pBRoute, currRoute });
 }

`RouteManager`则负责将`CurrentRoute`,`PersonalBestRoute`与`LocalBestRoute`各自的一部分组合形成新的`CurrentRoute`。

//updates a particle's velocity. The shorter the total distance the greater the velocity
public double UpdateVelocity(IRoute particleItinery, double weighting, double randomDouble)
{
    return (1 / particleItinery.TourDistance) * randomDouble * weighting;
}

//Selects a section of the route with a length  proportional to the particle's
// relative velocity.
public int GetSectionSize(IRoute particleItinery, double segmentVelocity, double totalVelocity)
{
    int length = particleItinery.DestinationIndex.Length;
    return Convert.ToInt32(Math.Floor((segmentVelocity / totalVelocity) * length));
}
  最后,`RouteManager`使用一个`RouteUpdater`来处理新的路径的构建。
//sets the selected BitArray mask so that
//only cities that have not been added already are available
//pointer is set to the start of the segment
public void SetSelectedMask(int pointer, IRoute section)
{
    int p = pointer;
    this.SelectedMask.SetAll(false);
    //foreach city in the section set the appropriate bit
    // in the selected mask
    for (int i = 0; i < section.SegmentSize; i++)
    {
        //set bit to signify that city is to be added if not already used
        this.SelectedMask[section.DestinationIndex[p]] = true;

        p++;
        //p is a circular pointer in that it moves from the end of the route
        // to the start
        if (p == section.DestinationIndex.Length)
        {
            p = 0;
        }
    }
    //in the AvailabilityMask, true=available, false= already used
    //remove cities from the SelectedMask  that have already been added
    this.SelectedMask.And(this.AvailabilityMask);
}

 //Updates the  new route by adding cities,sequentially from the route section
//providing the cities are not already present
public void SetDestinationIndex(int startPosition, IRoute section)
{


    int p = startPosition;
    for (int i = 0; i < section.SegmentSize; i++)
    {
        if (this.SelectedMask[section.DestinationIndex[p]])
        {
            this.destinationIndex[this.destinationIndexPointer] = section.DestinationIndex[p];
            this.destinationIndexPointer++;
        }
        p++;
        if (p == section.DestinationIndex.Length)
        {
            p = 0;
        }
    }
    //update the City AvailabilityMask
    //sets bits that represent cities that have been included to false
    this.AvailabilityMask.Xor(this.SelectedMask);
}

测试结果

(译者注:下面描述的参数和图中有一些出入,但原文就是这样写的。)

使用以下参数进行 100 个群体优化的测试:

1. 测试文件 Pr76DataSet.xml,包含 76 个城市,正确解决方案为 108,159 个。

2. 群体大小(粒子数量)为 80。

3. 每个群体优化的次数为 30000 次。

4. 组内“线人”数为 8 个 。

5. 重新分组之前的静态优化次数为 250 次。

6. 权重 W = 0.7, C1 = 1.4, C2 = 1.4。

测试结果:

1. 找到 7 个正解。 2. 最高错误率为 6%。 3. 平均误差为 2%。 4. 单个群的优化时间为 1 分 30 秒。

小结

粒子群优化算法可通过重复多次使用一个简单的算法来解决一些高度复杂的问题。

本文的版权归 StoneDemo 所有,如需转载请联系作者。

编辑于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏专知

【读书笔记】基于知识库的问答:生成查询图进行语义分析

【导读】将DBPedia和Freebase这样的大规模知识库组织并存储在一个结构化的数据库,这已成为支持开放领域问题问答的重要资源。 KB-QA的大多数方法基于...

4467
来自专栏AI研习社

MIT Taco 项目:自动生成张量计算的优化代码,深度学习加速效果提高 100 倍

我们生活在大数据的时代,但在实际应用中,大多数数据是 “稀疏的”。例如,如果用一个庞大的表格表示亚马逊所有客户与其所有产品的对应映射关系,购买某个产品以 “1”...

36111
来自专栏应用案例

关于机器学习,这可能是目前最全面最无痛的入门路径和资源!

之前搞机器学习的那帮人都喜欢用Python,所以Python慢慢就积攒了很多优秀的机器学习库,所谓的库,你就理解为别人封装好的一些具有某些功能的模块,我们可以通...

3077
来自专栏AI研习社

告别选择困难症,我来带你剖析这些深度学习框架基本原理

无论你喜欢或不喜欢,深度学习就在这里等着你来学习,伴随着技术淘金热而来的过多的可选项,让新手望而生畏。

1763
来自专栏生信技能树

如何通过Google来使用ggplot2可视化

今天是大年初二,这篇文章我只想传达一点: 没有什么菜鸟级别的生物信息学数据处理是不能通过Google得到解决方案的,如果有,请换个关键词继续Google! 第一...

3238
来自专栏AI科技评论

开发丨深度学习框架太抽象?其实不外乎这五大核心组件

许多初学者觉得深度学习框架抽象,虽然调用了几个函数/方法,计算了几个数学难题,但始终不能理解这些框架的全貌。 为了更好地认识深度学习框架,也为了给一些想要自己亲...

3694
来自专栏专知

【AlphaGo Zero 核心技术-深度强化学习教程代码实战06】给Agent添加记忆功能

【导读】Google DeepMind在Nature上发表最新论文,介绍了迄今最强最新的版本AlphaGo Zero,不使用人类先验知识,使用纯强化学习,将价值...

5146
来自专栏数据派THU

自然语言处理如何检查拼写错误?(Tensorflow实例教程、源代码)

原文:Towards Data Science 作者:Dave Currie 来源:机器人圈 本文长度为2400字,建议阅读5分钟 本文教你用TensorFlo...

1.1K8
来自专栏恰同学骚年

Unity3D游戏开发初探—2.初步了解3D模型基础

  简而言之,3D模型就是三维的、立体的模型,D是英文Dimensions的缩写。

1073
来自专栏专知

【AlphaGo Zero 核心技术-深度强化学习教程代码实战02】理解gym的建模思想

点击上方“专知”关注获取更多AI知识! 【导读】Google DeepMind在Nature上发表最新论文,介绍了迄今最强最新的版本AlphaGo Zero,不...

3415

扫码关注云+社区

领取腾讯云代金券