蚁群算法求包含34个国内城市的TSP，和最优解相差没那么大

文章目录

引言
蚁群觅食
算法原理
代码实现
- ACO求解TSP
- 整数规划求解TSP
相关阅读

引言

上一篇介绍的差分进化算法，很适合求解连续变量的优化问题；但针对组合优化问题，就不是很适用了。

至于哪一种智能优化算法更适合求解组合优化问题，我查阅了很多资料。不幸的是，虽然有很多公开测试数据集，但并没有令人信服的证据能表明某个智能优化算法是其中最好的。

即使如此，还是想找一个适合的智能优化算法来求解组合优化问题。一番权衡后，本文选择了蚁群算法，其中的原因包括以下两点：

蚁群算法的最初设计就是解决其中一大类组合优化问题，后续的应用几乎也都是组合优化问题，所以从天性上来说，两者是契合的；
蚁群算法最早是在1991年被提出的，而且持续有人在使用和研究，是一种经典的智能优化算法，所以具有学习的价值。

正文见下。

蚁群觅食

蚁群算法的设计借鉴了蚁群觅食的过程，因此为了能深刻地理解算法，有必要先了解清楚蚁群是如何觅食的。

研究发现，蚂蚁的食物搜索行为具有协作性。当蚂蚁在蚁穴和食物之间行走时，会释放一种被称为信息素的物质，这些物质形成一条指示轨迹。其他蚂蚁可以感知到环境中这种物质的存在和强度，在移动过程中，倾向于选择信息素强度高的轨迹。

以下图为例。 $t = 0$ 时刻，有20只蚂蚁位于A（蚁穴）处，要去往E（食物）处，另有20只蚂蚁位于E处，要返回A处。

$t = 1$ 时刻，20只蚂蚁从A到达B处的岔路口，因为此时没有信息素，所以10只蚂蚁走BCD路线，10只走BD路线，并分别在各自的路线上留下信息素；同理，从E出发的20只蚂蚁在D处也分成2个等数量的小组分别走DB和DCB路线，并留下信息素；另外，还有15只蚂蚁从A处开始出发。

到了 $t = 2$ 时，从E出发沿着DB路线的10只蚂蚁到达B处，而沿DCB行走的蚂蚁仅达到C处，因此此时BD路线上的信息素强度高于BCD路线；当另外15只从A出发的蚂蚁到达B处后，多数的蚂蚁（10只）会选择走BD路线，少量的蚂蚁（5只）选择走BCD路线。

这一正反馈机制会导致，BD路线上的信息素量相比BCD路线会越来越多，最终所有蚂蚁都找到BD这条最短路径。

算法原理

蚁群算法源于蚁群觅食，所以保留了一些类似的特征，比如个体间相互合作，通过信息素进行间接通信和正反馈机制等；但又高于蚁群觅食，主要体现在具有记忆能力，信息素的释放时机可以根据具体问题设定，允许局部优化、回退甚至预测能力等。

以下借助TSP问题（旅行商问题）实例，来详细描述蚁群算法的原理和和实现步骤。TSP问题指的是：旅行家要旅行n个城市，要求各个城市经历且仅经历一次，然后回到出发城市，并要求所走的路程最短。

（1）在蚁群算法的每一代循环中，每只蚂蚁都会独立地完成自己的TSP。蚂蚁完成TSP的过程中，不会选择已经访问过的城市，并且依赖状态转移概率选择下一个城市 $j$ 。状态转移概率的定义为：设第 $k$ 只蚂蚁当前所在的城市为 $i$ ，则从该城市移动到 $j$ 城市的状态转移概率为

其中， $\tau (i,j)$ 是城市对 $(i, j)$ 上的信息素，下一段再细说； $\eta (i,j)$ 是城市对 $(i, j)$ 上的启发信息，在TSP问题中， $\eta (i,j)$ 一般选取为城市间距离的倒数。 $N_k^i$ 是所有未经过的城市组成的集合， $\alpha,\beta$ 分别表示信息素和启发信息相对权重参数，它们控制 $\tau (i,j)$ 和 $\eta (i,j)$ 在决策中所占的比重。显然，信息素越多，距离越短，被选择的概率就越大。

（2）初始的信息素设置成所有城市对均相同即可。当所有蚂蚁完成了自己的TSP后，会根据自己选择的路径和路径总长度留下不同强度的信息素 $\tau (i,j)$ 。用 $\Delta\tau_{i,j}^k$ 表示第 $k$ 只蚂蚁在城市对 $(i, j)$ 上释放的量。在Ant Cycle System中，其值为

其中，$Q$是常数，$L_k$表示第$k$只蚂蚁在本次循环中所走路径的长度。显然，解的质量越好，在所经过的路径上留下的信息素$\Delta\tau_{i,j}^k$就越大。

所以，针对城市对 $(i, j)$ 上，信息素的总改变量为
$\Delta\tau_{i,j}=\sum_{k=1}^m \Delta\tau_{i,j}^k$
其中， $m$ 为蚂蚁数。

此外，算法中还引入了信息素会发机制。设信息素的保持系数为 $\rho$ ，则在完成一次完整循环后，信息素的值变为
$\tau_{i,j} = \rho\tau_{i,j} + \Delta\tau_{i,j}$

（3）判断是否达到退出条件，如果满足，则退出，否则便进入下一代循环。

代码实现

ACO求解TSP

求解随机生成的TSP实例，感觉有些太普通。所以决定加点难度，试着沿国内34个主要城市跑一圈，并把优化后的最短TSP可视化出来。

以下是Python代码。感恩网上大量的开源代码，让我完成可以少花很多时间来完成这些代码，同时保证TSP可视化后的颜值。

import math
import time

import numpy as np
import pandas as pd

from pyecharts import options as opts
from pyecharts.charts import Geo
from pyecharts.datasets import register_url
from pyecharts.globals import ChartType, SymbolType


def tsp_by_aco(D):
    # 城市数量
    city_cnt = D.shape[0]
    # 蚂蚁数量
    ant_cnt = 100
    # 迭代次数
    max_iter = 200
    # 信息素权重系数
    alpha = 1
    # 启发信息权重系数
    beta = 2
    # 信息素挥发速度
    rho = 0.1
    # 城市间球面距离
    distance = np.zeros((city_cnt, city_cnt))
    for i in range(city_cnt):
        for j in range(city_cnt):
            if i == j:
                # 相同城市不允许访问
                distance[i][j] = 1000000
            else:
                # 单位：km
                distance[i][j] = 6378.14 * math.acos(
                    math.cos(D[i][1]) * math.cos(D[j][1]) * math.cos(D[i][0] - D[j][0]) +
                    math.sin(D[i][1]) * math.sin(D[j][1]))
    # 启发信息，距离倒数
    eta = 1 / distance
    # 信息素矩阵
    tau = np.ones((city_cnt, city_cnt))
    # 每一代最优解
    generation_best_y = []
    generation_best_x = []
    # 种群
    population = np.zeros((ant_cnt, city_cnt)).astype(int)

    # 循环迭代
    for i in range(max_iter):
        # 城市转移概率
        prob_matrix = (tau ** alpha) * (eta ** beta)

        # TSP距离
        y = np.zeros((ant_cnt, 1))

        # 依次遍历每只蚂蚁
        for j in range(ant_cnt):
            # 设置TSP初始点为0
            population[j, 0] = 0

            # 选择后续城市
            for k in range(city_cnt - 1):
                # 已访问城市
                visit = set(population[j, :k + 1])
                # 未访问城市
                un_visit = list(set(range(city_cnt)) - visit)
                # 未访问城市转移概率归一化
                prob = prob_matrix[population[j, k], un_visit]
                prob = prob / prob.sum()
                # 轮盘赌策略选择下个城市
                next_point = np.random.choice(un_visit, size=1, p=prob)[0]
                # 添加被选择的城市
                population[j, k + 1] = next_point
                # 更新TSP距离
                y[j] += distance[population[j, k], population[j, k + 1]]
            # 更新TSP距离：最后一个城市->第0个城市
            y[j] += distance[population[j, -1], 0]

        # 保存当前代最优解
        best_index = y.argmin()
        generation_best_x.append(population[best_index, :])
        generation_best_y.append(y[best_index, :])

        # 计算信息素改变量，ACS模型，Q=1
        delta_tau = np.zeros((city_cnt, city_cnt))
        for j in range(ant_cnt):
            for k in range(city_cnt - 1):
                delta_tau[population[j, k], population[j, k + 1]] += 1 / y[j]
            delta_tau[population[j, city_cnt - 1], population[j, 0]] += 1 / y[j]

        # 信息素更新
        tau = (1 - rho) * tau + delta_tau

        print('iter: {}, best_f: {}'.format(i, generation_best_y[-1]))

    # 最优解位置
    best_generation_index = np.array(generation_best_y).argmin()

    return generation_best_x[best_generation_index], generation_best_y[best_generation_index]


def plot_tsp(cities):
    try:
        register_url("https://echarts-maps.github.io/echarts-countries-js/")
    except Exception:
        import ssl

        ssl._create_default_https_context = ssl._create_unverified_context
        register_url("https://echarts-maps.github.io/echarts-countries-js/")
    '''
    https://echarts-maps.github.io/echarts-countries-js/preview.html
    这个网站上显示的各个国家中文名称， 可以写在下面的maptype里面
    '''
    title1 = "中国主要城市TSP_by_ACO"

    c = (
        Geo().add_schema(
            maptype='china',  # 可以输入国家名字，比如"瑞士"
            itemstyle_opts=opts.ItemStyleOpts(color='#323c48', border_color='#111'), )  # 设置地图颜色和边框色
    )

    for i in range(len(cities)):
        c.add(
            "城市",  # 第一个add数据的标题
            [(cities['城市'].iloc[i], i)],
            type_=ChartType.EFFECT_SCATTER,  # 使用点的样式，并设置点的颜色，点的大小都是一样的！
            symbol_size=6,  # 设置点的大小
            color='red', )  # 点的颜色

    for i in range(len(cities) - 1):
        c.add(
            "tsp路线",
            [(cities['城市'].iloc[i], cities['城市'].iloc[i + 1])],  # 城市顺序
            type_=ChartType.LINES,
            effect_opts=opts.EffectOpts(
                symbol=SymbolType.ARROW, symbol_size=6, color='yellow'),  # 线上的小箭头的颜色
            linestyle_opts=opts.LineStyleOpts(curve=0.2))  # 设置两点间线缆的弯曲度
        c.set_series_opts(label_opts=opts.LabelOpts(is_show=False))
        c.set_global_opts(title_opts=opts.TitleOpts(title=title1),
                          toolbox_opts=opts.ToolboxOpts())

    c.render('tsp_solution_by_aco.html')


if __name__ == '__main__':
    # 城市及其经纬度
    original_cities = [['西宁', 101.74, 36.56],
                       ['兰州', 103.73, 36.03],
                       ['银川', 106.27, 38.47],
                       ['西安', 108.95, 34.27],
                       ['郑州', 113.65, 34.76],
                       ['济南', 117, 36.65],
                       ['石家庄', 114.48, 38.03],
                       ['太原', 112.53, 37.87],
                       ['呼和浩特', 111.65, 40.82],
                       ['北京', 116.407526, 39.90403],
                       ['天津', 117.200983, 39.084158],
                       ['沈阳', 123.38, 41.8],
                       ['长春', 125.35, 43.88],
                       ['哈尔滨', 126.63, 45.75],
                       ['上海', 121.473701, 31.230416],
                       ['杭州', 120.19, 30.26],
                       ['南京', 118.78, 32.04],
                       ['合肥', 117.27, 31.86],
                       ['武汉', 114.31, 30.52],
                       ['长沙', 113, 28.21],
                       ['南昌', 115.89, 28.68],
                       ['福州', 119.3, 26.08],
                       ['台北', 121.3, 25.03],
                       ['香港', 114.173355, 22.320048],
                       ['澳门', 113.54909, 22.198951],
                       ['广州', 113.23, 23.16],
                       ['海口', 110.35, 20.02],
                       ['南宁', 108.33, 22.84],
                       ['贵阳', 106.71, 26.57],
                       ['重庆', 106.551556, 29.563009],
                       ['成都', 104.06, 30.67],
                       ['昆明', 102.73, 25.04],
                       ['拉萨', 91.11, 29.97],
                       ['乌鲁木齐', 87.68, 43.77]]
    original_cities = pd.DataFrame(original_cities, columns=['城市', '经度', '纬度'])
    # 使用ACO算法求解TSP
    time0 = time.time()
    best_x, best_y = tsp_by_aco(original_cities[['经度', '纬度']].values * math.pi / 180)
    print('使用ACO求解TSP，耗时: {} s'.format(time.time() - time0))
    # 按最优解顺序，生成访问城市次序
    original_cities['seq'] = best_x
    sort_cities = original_cities.sort_values(by='seq')
    sort_cities.loc[len(sort_cities)] = sort_cities.loc[0]
    # 绘制TSP路径
    plot_tsp(sort_cities)

运行代码后，能够得到一个网页版的动图，可以表示运行的轨迹，这里截取了一帧如下。算法输出的最短路径是15944.43815449，从图里其实也能看出来，这个解应该不是全局最优解。需要说明的是，每次运行的结果大概率是不一样的，所以如果大家复跑了这段代码，最优解的值和图的形状不一样是正常的。

整数规划求解TSP

之前的文章里介绍过整数规划可以求解全局最优解，所以本节再用这个方法解一遍TSP，然后把全局最优解拿来和ACO算法的结果做一下对比。TSP问题的整数规划建模不是本文重点，这里就不多描述了，有兴趣的可以参考：优化│TSP中两种不同消除子环路的方法及callback实现（Python调用Gurobi求解）。这里直接给出使用ortools求解模型的代码。

import pandas as pd
import math
import numpy as np
import time
from ortools.linear_solver import pywraplp

from pyecharts import options as opts
from pyecharts.charts import Geo
from pyecharts.datasets import register_url
from pyecharts.globals import ChartType, SymbolType


def calc_by_ortools(data):
    # 生成ortools求解器，使用SCIP算法
    solver = pywraplp.Solver.CreateSolver('SCIP')

    # 优化变量
    x = {}
    u = {}
    for i in range(len(data)):
        u[i] = solver.NumVar(0, solver.infinity(), 'u[%i]' % i)
        for j in range(len(data)):
            x[i, j] = solver.IntVar(0, 1, 'x[%i,%i]' % (i, j))

    # 目标函数:
    obj_expr = [data[i][j] * x[i, j] for i in range(len(data)) for j in range(len(data))]
    solver.Minimize(solver.Sum(obj_expr))

    # 约束条件一: 出度为1
    for i in range(len(data)):
        constraint_expr = [x[i, j] for j in range(len(data))]
        solver.Add(solver.Sum(constraint_expr) == 1)
    # 约束条件二: 入度为1
    for j in range(len(data)):
        constraint_expr = [x[i, j] for i in range(len(data))]
        solver.Add(solver.Sum(constraint_expr) == 1)
    # 约束条件三: 消除子环
    for i in range(1, len(data)):
        for j in range(1, len(data)):
            if i != j:
                solver.Add(u[i] - u[j] + len(data) * x[i, j] <= len(data) - 1)

    # 模型求解
    status = solver.Solve()

    # 模型求解成功, 打印结果
    if status == pywraplp.Solver.OPTIMAL:

        # 最优目标函数值
        print('best_f =', solver.Objective().Value())

        # 最优次序
        cnt = 0
        print(0, end=' ')
        k = 0
        best_x = [0]
        while cnt < len(data):
            for j in range(len(data)):
                if x[k, j].solution_value() == 0:
                    continue
                print('->', j, end=' ')
                k = j
                cnt += 1
                best_x.append(j)
                break
        print('\n')
        return best_x

    else:
        print('not converge.')
        return []


def plot_optimal_routes(cities):
    try:
        register_url("https://echarts-maps.github.io/echarts-countries-js/")
    except Exception:
        import ssl

        ssl._create_default_https_context = ssl._create_unverified_context
        register_url("https://echarts-maps.github.io/echarts-countries-js/")
    '''
    https://echarts-maps.github.io/echarts-countries-js/preview.html
    这个网站上显示的各个国家中文名称， 可以写在下面的maptype里面
    '''
    title1 = "中国主要城市TSP_by_ortools"

    c = (
        Geo()
            .add_schema(
            maptype='china',  # 可以输入国家名字，比如"瑞士"
            itemstyle_opts=opts.ItemStyleOpts(color='#323c48', border_color='#111'), )  # 设置地图颜色和边框色
    )

    for i in range(len(cities)):
        c.add(
            "城市",  # 第一个add数据的标题
            [(cities['城市'].iloc[i], i)],
            type_=ChartType.EFFECT_SCATTER,  # 使用点的样式，并设置点的颜色，点的大小都是一样的！
            symbol_size=6,  # 设置点的大小
            color='red', )  # 点的颜色

    for i in range(len(cities) - 1):
        c.add(
            "tsp路线",
            [(cities['城市'].iloc[i], cities['城市'].iloc[i + 1])],  # 城市顺序
            type_=ChartType.LINES,
            effect_opts=opts.EffectOpts(
                symbol=SymbolType.ARROW, symbol_size=6, color='yellow'),  # 线上的小箭头的颜色
            linestyle_opts=opts.LineStyleOpts(curve=0.2))  # 设置两点间线缆的弯曲度
        c.set_series_opts(label_opts=opts.LabelOpts(is_show=False))
        c.set_global_opts(title_opts=opts.TitleOpts(title=title1),
                          toolbox_opts=opts.ToolboxOpts())

    c.render('tsp_solution_by_ortools.html')


if __name__ == '__main__':
    # 城市及其经纬度
    original_cities = [['西宁', 101.74, 36.56],
                       ['兰州', 103.73, 36.03],
                       ['银川', 106.27, 38.47],
                       ['西安', 108.95, 34.27],
                       ['郑州', 113.65, 34.76],
                       ['济南', 117, 36.65],
                       ['石家庄', 114.48, 38.03],
                       ['太原', 112.53, 37.87],
                       ['呼和浩特', 111.65, 40.82],
                       ['北京', 116.407526, 39.90403],
                       ['天津', 117.200983, 39.084158],
                       ['沈阳', 123.38, 41.8],
                       ['长春', 125.35, 43.88],
                       ['哈尔滨', 126.63, 45.75],
                       ['上海', 121.473701, 31.230416],
                       ['杭州', 120.19, 30.26],
                       ['南京', 118.78, 32.04],
                       ['合肥', 117.27, 31.86],
                       ['武汉', 114.31, 30.52],
                       ['长沙', 113, 28.21],
                       ['南昌', 115.89, 28.68],
                       ['福州', 119.3, 26.08],
                       ['台北', 121.3, 25.03],
                       ['香港', 114.173355, 22.320048],
                       ['澳门', 113.54909, 22.198951],
                       ['广州', 113.23, 23.16],
                       ['海口', 110.35, 20.02],
                       ['南宁', 108.33, 22.84],
                       ['贵阳', 106.71, 26.57],
                       ['重庆', 106.551556, 29.563009],
                       ['成都', 104.06, 30.67],
                       ['昆明', 102.73, 25.04],
                       ['拉萨', 91.11, 29.97],
                       ['乌鲁木齐', 87.68, 43.77]]
    original_cities = pd.DataFrame(original_cities, columns=['城市', '经度', '纬度'])
    D = original_cities[['经度', '纬度']].values * math.pi / 180

    # 城市间球面距离
    city_distance = np.zeros((len(original_cities), len(original_cities)))
    for i in range(len(original_cities)):
        for j in range(len(original_cities)):
            if i == j:
                city_distance[i][j] = 100000
            else:
                city_distance[i][j] = 6378.14 * math.acos(
                    math.cos(D[i][1]) * math.cos(D[j][1]) * math.cos(D[i][0] - D[j][0]) +
                    math.sin(D[i][1]) * math.sin(D[j][1]))

    # 使用整数规划求解TSP
    time0 = time.time()
    best_x = calc_by_ortools(city_distance)
    print('使用ortools求解TSP，耗时: {} s'.format(time.time() - time0))

    # 绘制TSP路径
    if len(best_x) > 0:
        best_routes = []
        for i in range(len(best_x)):
            best_routes.append(original_cities['城市'].iloc[best_x[i]])
        best_routes = pd.DataFrame(best_routes, columns=['城市'])
        plot_optimal_routes(best_routes)

运行程序后，得到最优解为15614.849985682848，可视化路径如下图所示，这个图里的路径看上去更符合我们对全局最优解的认知。相比之下，蚁群算法的结果距离最优解的GAP差不多是2%，看起来也没差那么多。

best_f = 15614.849985682848
0 -> 1 -> 2 -> 3 -> 4 -> 6 -> 7 -> 8 -> 13 -> 12 -> 11 -> 9 -> 10 -> 5 -> 17 -> 16 -> 14 -> 15 -> 22 -> 21 -> 20 -> 18 -> 19 -> 25 -> 23 -> 24 -> 26 -> 27 -> 31 -> 28 -> 29 -> 30 -> 32 -> 33 -> 0 
使用ortools求解TSP，耗时: 88.28608894348145 s