LK算法、LKH算法介绍及Python实现

原创 2023-11-23 20:55:16 1.1w 阅读 AI写同款

收录于

优化算法

QQVQQ_

码龄8年墨衍会员 · AI 创作全网分发

2610粉丝107篇原创

关注

LKH求解TSP 文章目录一：预备知识LK算法求解TSP二：介绍三：改进的LK算法（即LKH1）一：预备知识 LK算法求解TSP 二：介绍 LK算法通常被考虑是一个有效的办法来对对称的TSP问题产生一个最优或近似最优解。然后，LK算法的执行并不简单。LK算法中有需要操作和步骤需要完成，且这些步骤在性能上都起到了很大的作用。 LKH算法（改进的LK算法）和LK算法有很多细节上的不同，其中最明显的不同在于搜索策略的不同。LKH算法使用比LK算法更大更复杂的搜索步骤，并且使用sensitivity analysis来指导和限制搜阅读详情

AI摘要

本文详细解析了LK算法的基本步骤，相关知识点，以及正增益、不相交等准则。介绍了LKH算法的改进规则，并提供了Python代码实现，探讨了如何通过启发式规则加速搜索。适合理解和实践TSP问题的旅行商算法。

0. 前言

有关算法的详细信息参考下面这两个网址：

Extreme Algorithmshttp://www.seas.gwu.edu/~simhaweb/champalg/tsp/tsp.html

LKH (Keld Helsgaun)http://www.akira.ruc.dk/~keld/research/LKH/

1. 基本LK算法

LK算法是基于交换转换的实质性概括。假设T是一个非最优但可行的解决方案。那么我们可以说T不是最优的，因为T中有k个元素 $x_{1},...,x_{k}$ 是“不合适”；为了使T最优，它们应该被 S-T 的k元素 $y_{1},...,y_{k}$ 代替，问题只是识别 k 和 x 和 y 。

由于我们不知道k应该是什么，所以修复k然后考虑这个固定k的所有可能的T的k-subsets似乎是人为的。相反，我们尝试逐个元素找到 k 和 $x_{1},...,x_{k}$ 和 $y_{1},...,y_{k}$ 作为最好的。因此，我们将首先尝试将 $x_{1}$ 和 $y_{1}$ 识别为“最不合适”的对；然后，选择 $x_{1}$ 和 $y_{1}$ 并暂时搁置，我们发现 $x_{2}$ 和 $y_{2}$ 是其余集合中最不合适的一对；等等。

1.1 基本步骤

1.2 相关知识点

设 $\left | x_{i} \right |,\left | x_{j} \right |$ 表示 $x_{i}$ 和 $y_{i}$ 的长度，定义 $g_{i} = \left | x_{i} \right | - \left | y_{i} \right |$

LK算法基于一个重要的simple fact：如果一系列序列具有正和，则这些序列存在循环排列，使得每个部分和都是正的。

证明：

设 k 是使得 $g_{1} + ... + g_{k-1}$ 最小的索引，则

特别是，由于我们正在寻找具有正和的 $g_{i}$ 序列，我们只需要考虑部分和总是正的增益序列。这种增益标准使我们能够大大减少我们需要检查的序列数量；这是算法停止规则的核心。

1.3 算法的伪代码

算法的有关细节如下：

Step1. 随机生成一个初始旅行 $T$ ，即随机生成一个可行解。

Step2. $i$ 用于记录交换边的数量， $t_{1}$ 为在当前 $T$ 状态下，未选择的一个城市节点。

Step3. 在 $T$ 中选择一条边 $x_{1}=(t_{1},t_{2})$ 。当 $t_{1}$ 确定时， $t_{2}$ 有两个选择,即 $x_{1}$ 也仅有两个选择。两个选择都应该进行判断，寻找潜在的最优解。

Step4. 当 $x_{1}=(t_{1},t_{2})$ 确定时，则 $y_{1}=(t_{2},t_{3})$ 的一个端点确定，此时 $t_{3}$ 不应该等于 $t_{2}$ ，此时选择的 $y_{1}$ 应当保证 $G_{1}>0$ , 如果没有满足条件的应当尝试下一个 $x_{1}$ ，如果 $x_{1}$ 为当前 $t_{1}$ 下最后一个满足条件的 $x_{1}$ 则尝试下一个 $t_{1}$ 。

Step5. 令 $i= i+1$ ，进行下一组边的尝试。

Step6. 对于给定的 $y_{i-1}$ , $x_{i}$ 有两种选择，但是两种选择中仅有一种可以闭合整个旅行，其中一个无法闭合，仅当选择 $i = 2$ ，即选择 $x_{2}$ 时，才允许不闭合整个旅行的存在，这可以增大获得全局最优解的可能性。 $x_{i}$ ≠ $y_{s}$ , $s<i$ . 意味着 $x_{i}$ 不能是先前已经断开的边。

Step7. 在 $x_{i}=(t_{2i-1},t_{2i})$ 确定时，则 $y_{i}=(t_{2i},t_{2i+1})$ 的 $t_{2i}$ 已经确定，此时需要确定 $t_{2i+1}$ ,且此时 $y_{i}$ 的存在必须确保 $x_{i+1}$ 也要存在， $x_{i+1}$ 存在的条件参考Step6中的（a）、（b）。如果没有满足条件的则更换 $y_{2}$ 进行尝试，如果有满足条件的则继续增加交换的边，看是否可以继续增加。

Step8. 步骤8是对所有 $y_{2}$ 的遍历

Step9. 步骤9是对所有 $x_{2}$ 的遍历

Step10. 步骤10是对所有 $y_{1}$ 的遍历

Step11. 步骤11是对所有 $x_{1}$ 的遍历

Step12. 步骤12是对所有 $t_{1}$ 的遍历

Step13. 停止条件。

1.4 算法需要满足的准则

1. 准则一（顺序交换准则）

$y_{i}$ 与 $x_{i}$ 、 $x_{i+1}$ 分别共享一个节点。如果 $t_{1}$ 是 $x_{1}$ 其中一个节点，则当 $i\ge 1$ 时， $\mathrm{x}_{\mathrm{i}}=\left(\mathrm{t}_{2 \mathrm{i}-1}, \mathrm{t}_{2 \mathrm{i}}\right), \mathrm{y}_{\mathrm{i}}=\left(\mathrm{t}_{2 \mathrm{i}}, \mathrm{t}_{2 \mathrm{i}+1}\right)$ ， $\mathrm{x}_{\mathrm{i}+1}=\left(\mathrm{t}_{2 \mathrm{i}+1}, \mathrm{t}_{2 \mathrm{i}+2}\right)$ ，如下图所示

如图， $\left(x_{1}, y_{1}, x_{2}, y_{2}, x_{3}, \ldots, x_{r}, y_{r}\right)$ 构成相邻连接链。

集合X与集合Y交换产生旅行的必要（但不充分）条件是连接链被关闭，即 $\mathrm{y}_{\mathrm{r}}=\left(\mathrm{t}_{2 \mathrm{r}}, \mathrm{t}_{1}\right)$ 。这种交换称为顺序交换。

通常，可以通过适当编号受影响链路作为顺序交换来实现旅行的改进。但是，情况并非总是如此。下图显示了一个无法进行顺序交换的示例。

2. 准则二（可行性准则）

对于选择的 ${x}_{\mathrm{i}}=\left(\mathrm{t}_{2 \mathrm{i}-1}, \mathrm{t}_{2 \mathrm{i}}\right)$ ,当 $t_{2i}$ 连接到 $t_{1}$ 时，确保得到的结果是一个可行的旅行，这个准则用于 ${i} \geq 3$ 确保可以关闭旅行，形成一个闭环。

3. 准则三（正增益准则）

要求选择的 $y_{i}$ 始终可以确保所提出的交换集合的增益 $G_{i}$ 为正，假设 $\mathrm{g}_{\mathrm{i}}=\mathrm{c}\left(\mathrm{x}_{\mathrm{i}}\right)-\mathrm{c}\left(\mathrm{y}_{\mathrm{i}}\right)$ 是 $x_{i}$ 与 $y_{i}$ 的收益，则 $\mathrm{G}_{\mathrm{i}} = \mathrm{g}_{1}+\mathrm{g}_{2}+\ldots+\mathrm{g}_{\mathrm{i}}$ 。该准则在算法的效率中起着重要作用。

4. 准则四（不相交准则）

最后，要求集合X和Y是不相交的。这简化了编码，也减少了运行时间并提供了有效的停止标准。

1.5 基本LK算法的改进

为了限制搜索，Lin和Kernighan通过引入以下规则来改进算法：

5. 准则五

在寻找一个进入 $T$ 的边时， $y_{i}=\left(t_{2 i}, t_{2 i+1}\right)$ 中的 $t_{2i}$ 被限制为最近的5个邻居。

6. 准则六

对于i≥4，如果 $T$ 中的任何一个边 $x_{i}$ 是少数（2-5个）可行解的共同环节，则该环节必须被打破。

7. 准则七

如果当前旅游与之前的解决方案旅游相同，则停止搜索。

规则5和6是启发式规则。它们是基于对哪些链接可能属于最优旅游的预期。它们可以节省运行时间，但有时是以不能获得最佳解决方案为代价的。

规则7也可以节省运行时间，但对找到的解决方案的质量没有影响。如果一个游程与之前的解决方案游程相同，那么试图进一步改进它就没有意义了。因此，检查是否有更多改进可能所需的时间（检查时间）可以被节省。根据Lin和Kernighan的说法，以这种方式节省的时间通常是运行时间的30%到50%。

除了这些主要是为了限制搜索的细化之外，Lin和Kernighan还增加了一些细化，其目的主要是为了指导搜索。在算法可以选择备选方案的情况下，启发式规则被用来给这些备选方案赋予优先权。在只必须选择其中一个替代方案的情况下，选择具有最高优先权的方案。在必须尝试几个备选方案的情况下，这些备选方案按优先级递减的顺序被尝试（使用回溯法）。更具体地说，使用了以下规则。

8. 准则八

当边 $y_{i}(i \geq 2)$ 被选择时，每个可能的选择都被赋予优先级 $c\left(x_{i+1}\right)-c\left(y_{i}\right)$ 。

9. 准则九

如果 $x_{4}$ 有两种选择，则选择 $c\left(x_{4}\right)$ 最高的一种。

2. LKH算法

待完成

3. Python实现

此实现仅仅是一个简单版本，还需要进一步优化代码结构

from random import choice
import math
import matplotlib.pyplot as plt


def get_distances(coordinates):
    """根据城市坐标获取城市距离，返回关于城市距离的矩阵"""
    city_distance = []
    city_number_ = len(coordinates)
    for i in range(city_number_):
        city_i = []
        i_x, i_y = coordinates[i][0], coordinates[i][1]
        for j in range(city_number_):
            j_x, j_y = coordinates[j][0], coordinates[j][1]
            city_i.append(math.sqrt((i_x - j_x) ** 2 + (i_y - j_y) ** 2))
        city_distance.append(city_i)
    return city_distance


def get_matrix(city_number_, file_path, row):
    """读取.tsp文件，并返回距离矩阵"""
    f = open(file_path, "r")  # 以只读模式打开文件
    lines = f.readlines()  # 获取所有行并保存到lines列表中
    datas = []  # 保存数据
    for line in lines[row - 1:-1]:  # 读取指定的所有行
        line.strip('\n')  # 去掉换行符
        datas.extend(line.split())  # 去掉空格
    f.close()  # 关闭文件

    dis_matrix = []
    # 先获取下三角矩阵
    line = []
    for data in datas:
        value = int(data)
        line.append(value)
        if value == 0:
            dis_matrix.append(line[:])
            line = []
    # 再对下三角矩阵进行修改
    for i in range(city_number_):
        for j in range(i + 1, city_number_):
            dis_matrix[i].append(dis_matrix[j][i])
    return dis_matrix


def initialize(number):
    """生成初始解，number为城市数量"""
    # city_list = [v for v in range(number)]
    # path = []
    # for i in range(number):
    #     city = choice(city_list)
    #     city_list.remove(city)
    #     path.append(city)
    path = list(range(number))
    return path


def get_path_fitness(path):
    """计算路径path的距离"""
    path_distance = distance_matrix[path[0]][path[-1]]
    for i in range(len(path) - 1):
        path_distance += distance_matrix[path[i]][path[i + 1]]
    return path_distance


def run_algorithm(N):
    """执行算法，N为算法运行次数"""
    iterate_result = []  # 每次迭代的结果
    best_value = 99999
    best_path = []
    for n in range(N):
        print(f"========================执行第{n}次==========================")
        path = LK()
        fitness = get_path_fitness(path)
        iterate_result.append(fitness)
        if fitness < best_value:
            best_path = path
            best_value = fitness
    return best_value, best_path, iterate_result


def transform_path(path):
    """将基于位置表示的路径转化为基于边的表示的路径"""
    T = []
    for i in range(len(path) - 1):
        T.append((path[i], path[i + 1]))
    T.append((path[-1], path[0]))
    return T


def transform_T(T):
    """将基于边表示的路径转化为基于位置表示的路径"""
    path = [T[0][0], T[0][1]]
    length = len(T)
    del T[0]
    for i in range(length - 2):
        for t in T:
            if t[0] == path[-1]:
                path.append(t[1])
                T.remove(t)
                break
            elif t[1] == path[-1]:
                path.append(t[0])
                T.remove(t)
                break
    return path


def LK():
    """LK算法"""
    p = False

    # 初始化路径
    path = initialize(city_number)
    T = transform_path(path)
    n = 0

    # 第一层循环，遍历所有的t1（注意，当T更改时需要t1要重新开始循环）
    while n < city_number:
        t1 = n
        x1_list = get_x1_list(t1, T)

        # 第二层循环，遍历t1对应的所有x1，实际x1最多两种可能性
        for x1 in x1_list:
            t2 = (set(x1) - {t1}).pop()

            # 第三层循环，遍历x1对应的所有y1
            y1_list = get_y1_list(t2, x1, T)
            if not y1_list:  # 没有满足条件的y1,进行下一个x1判断
                print(f"x1={x1}时没有满足条件的y1 ")
                continue
            for y1 in y1_list:

                # 第四层循环，遍历y1对应的所有x2,其中x2最多两种可能性，其中一种选择无法构成回路
                p, x2_list = get_x2_list(x1, y1, T)
                if p:  # 已存在较优的回路，更改T重新选择t1
                    break
                for x2 in x2_list:

                    # 第五层循环，遍历x2对应的所有y2
                    y2_list = get_y2_list(y1, [x1, x2], T)
                    if y2_list:  # 存在使增益增加的y2,且x3存在
                        for y2 in y2_list:

                            # 进入r_opt搜索
                            print(f"===========================开始r_opt=======================")
                            print(f"||t1={t1}||x1={x1}||y1={y1}||x2={x2}||y2={y2}")
                            p, p1, T = r_opt(T, [x1, x2], [y1, y2])
                            print("===========================结束r_opt=======================")
                            if p:
                                break
                            if p1:  # 进入下一个y2的尝试
                                continue
                    else:  # 不存在满足条件的y2
                        continue  # 直接进入下一循环，重新选择x2
                    if p:
                        break
                if p:
                    break
            if p:
                n = 0
                break
        if p:
            n = -1  # T被重置，因此需要t1也要重新开始
            p = False
            print("存在改进，重置t1起点")
        n += 1
    return transform_T(T)


def r_opt(T, tried_x, tried_y):
    """对T进行r_opt优化，tried_x=[x1,x2]，tried_y=[y1,y2]"""
    # 当y2=（t4,t5）确定时，x3则只有两种选择包含t4或者t5
    # 寻找所有的x3
    p = False  # 寻找下一个t1的标志
    p1 = False  # 寻找下一个y2的标志
    i = 2
    while True:
        i += 1
        print(f"进行第{i}对边搜索")
        # 选择xi
        p6, xi, T = choose_xi(T, tried_x, tried_y)  # 进入此步意味着x(i+1)必然存在，但不一定使T更好
        if p6:  # T可以被改善,且在choose_xi()函数中已经对T进行了修改
            p = True
            break

        # 选择yi
        tried_x.append(xi)
        p7, yi = choose_yi(T, tried_x, tried_y)
        if not p7:  # 没有满足条件的yi,尝试下一个y2
            p1 = True
            break
        tried_y.append(yi)
    return p, p1, T


def choose_xi(T, tried_x, tried_y):
    """如果选择的t_2i连接t1后,构成新回路的收益大于旧回路收益则更改旧回路返回True"""
    # yi和xi+1必须共享一个端点，即yi_1与xi也必须共享一个端点
    # yi=（t2i，t2i+1）和xi+1=（t2i+1，t2i+2）

    t1 = (set(tried_x[0]) - set(tried_y[0])).pop()
    t2i_1 = list(set(tried_y[-1]) - set(tried_x[-1]))[0]
    xi = (0, 0)
    p = False
    for t in T:  # xi为二选一
        t2i = (set(t) - {t2i_1}).pop()
        if (t2i_1 in t) and (t not in tried_y) and (t not in tried_x) and (t2i != t1):
            yi = (t1, t2i)
            X = tried_x[:]
            X.append(t)
            Y = tried_y[:]
            Y.append(yi)
            if judge_swap(T, X, Y):  # 可以进行边的交换
                xi = t
                if judge_fitness(X, Y):
                    T = list((set(T) - set(X)) | set(Y))
                    print(f"xi存在改进，改进后的T为\n{T}")
                    p = True
                break

    return p, xi, T


def choose_yi(T, tried_x, tried_y):
    """尝试寻找满足条件的yi，如果找到则返回yi和True"""
    # xi = (t2i-1,t2i)
    # yi = (t2i,t2i+1)
    # 当xi确定时，yi的某个端点则确定

    t1 = (set(tried_x[0]) - set(tried_y[0])).pop()
    t2i = (set(tried_x[-1]) - set(tried_y[-1])).pop()

    for city in range(city_number):  # 遍历所有城市寻找满足条件的yi
        yi = (t2i, city)
        X = tried_x[:]
        Y = tried_y[:]
        Y.append(yi)
        G_i = judge_fitness(X, Y)  # 收益值
        if (G_i > 0) and (yi not in T) and ((city, t2i) not in T) and (city != t2i) and (city != t1):
            # 接下来判断xi+1是否存在，xi+1=(city，t_)
            for x in T:
                if city in x:
                    t_ = (set(x) - {city}).pop()
                    X = tried_x[:]
                    X.append(x)
                    Y = tried_y[:]
                    Y.extend([yi, (t_, t1)])
                    if judge_swap(T, X, Y) and (x not in tried_x):
                        return True, yi

    return False, (0, 0)


def judge_swap(T, X, Y):
    """判断交换能能否构成闭环，可以则返回True"""
    T_ = (set(T) - set(X)) | set(Y)
    t = 0
    i = 0
    while i <= len(T):
        for x in T_:
            if t in x:
                i += 1
                x_list = list(x)
                x_list.remove(t)
                t = x_list[0]
                if (t == 0) and (i < city_number - 1):
                    return False
    return True


def judge_fitness(X, Y):
    p = False
    fitness = 0
    for i in range(len(X)):
        fitness += distance_matrix[X[i][0]][X[i][1]] - distance_matrix[Y[i][0]][Y[i][1]]
    if fitness > 0:
        p = True
    return p


def get_T_fitness(T):
    """获取T的路径距离"""
    distance = 0
    for t in T:
        distance += distance_matrix[t[0]] + distance_matrix[t[1]]
    return distance


def get_x1_list(t1, T):
    """获取t1对应的所有可能的x1"""
    x1_list = []
    for v in T:
        if t1 in v:
            x1_list.append(v)
    return x1_list


def get_y1_list(t2, x1, T):
    """获取在x1确定下的所有的y1"""
    y1_list = []
    for city in range(city_number):
        G1 = distance_matrix[x1[0]][x1[1]] - distance_matrix[t2][city]
        if ((t2, city) not in T) and ((city, t2) not in T) and (G1 > 0) and (t2 != city):
            y1_list.append((t2, city))
    return y1_list


def get_x2_list(x1, y1, T):
    """x2有两种选择，只有一种可以形成回路"""
    t1 = (set(x1) - set(y1)).pop()
    p = False
    x2_list = []

    t2 = (set(x1) - {t1}).pop()
    t3 = (set(y1) - {t2}).pop()
    for x in T:
        if (t3 in x) and (t1 not in x):
            x2_list.append(x)
            t4 = (set(x) - {t3}).pop()
            X = [x1, x]
            Y = [y1, (t4, t1)]
            if judge_swap(T, X, Y) and (judge_fitness(X, Y)):
                print("！！！")
                T.remove(x)
                T.remove(x1)
                T.extend(Y)
                print(f"在x2处可改进，改进后的T为：\n{T}")
                p = True
                break
    return p, x2_list


def get_y2_list(y1, tried_x, T):
    """尝试寻找满足条件的y2，如果找到则返回所有y2列表"""
    t1 = (set(tried_x[0]) - set(y1)).pop()
    t4 = (set(tried_x[1]) - set(y1)).pop()

    y2_list = []
    # 遍历所有城市寻找满足条件的y2
    for city in range(city_number):
        y2 = (t4, city)
        G_i = judge_fitness(tried_x, [y1, y2])  # 收益值
        if (G_i > 0) and (y2 not in T) and ((city, t4) not in T) and (t4 != city) and (city != t1):
            # 接下来判断x3是否存在，x3=(city，t6)
            for x in T:
                if city in x:
                    t6 = (set(x) - {city}).pop()
                    X = tried_x[:]
                    X.append(x)
                    Y = [y1, y2, (t6, t1)]
                    if judge_swap(T, X, Y):
                        y2_list.append(y2)
    return y2_list


def get_neighbors(n):
    """获取所有城市的前n个邻居"""
    neighbors = {}
    for i in range(len(distance_matrix)):
        temp_list = sorted(distance_matrix[i])
        neighbors[i] = [distance_matrix[i].index(temp_list[v + 1]) for v in range(n)]
    return neighbors


if __name__ == '__main__':
    # 城市距离矩阵
    # distance_matrix = [
    #     [0, 64, 378, 519, 434, 200],
    #     [64, 0, 270, 455, 375, 164],
    #     [378, 270, 0, 70, 265, 344],
    #     [519, 455, 70, 0, 223, 428],
    #     [434, 375, 265, 223, 0, 173],
    #     [200, 164, 344, 428, 173, 0]
    # ]  # 6个城市，最优1248
    distance_matrix = get_matrix(24, "gr24.tsp", 8)
    city_number = len(distance_matrix)
    Run_Number = 1  # 迭代次数

    print("开始执行算法！！！")
    result = list(run_algorithm(Run_Number))  # 执行算法
    print("=======================================================")
    print(f"最优值为：{result[0]}，对应路径为：《{result[1]}》")

    # a = judge_swap([(6, 7), (9, 1), (11, 8), (4, 14), (13, 5), (13, 7), (15, 8), (14, 0), (9, 4), (12, 16), (1, 10),
    #                 (5, 2), (6, 16), (12, 3), (10, 2), (11, 3), (15, 0)], [(6, 7)], [(6, 7)])
    # print(a)

目前仅实现基本的LK算法

接下来将进行LKH算法的实现

参考文献：

[1] Helsgaun K . An effective implementation of the Lin–Kernighan traveling salesman heuristic[J]. European Journal of Operational Research, 2000, 126(1):106-130.

国内网站有关该算法介绍的文章见下面的网址：

TSP 问题的 LKH 方法介绍_xlh9718的博客-优快云博客_lkh算法

LKH求解TSP_yang_rookie的博客-优快云博客_lkh算法

请问有谁了解旅行商问题的LKH算法？可以交流一下吗？ - 知乎

tsp的理论和实践(9)LK简介 - 掘金

标签

#python #算法 #启发式算法 #图搜索算法