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

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

目录

0. 前言

1. 基本LK算法

1.1 基本步骤

1.2 相关知识点

1.3 算法的伪代码

1.4 算法需要满足的准则

1. 准则一(顺序交换准则)

2. 准则二(可行性准则)

3. 准则三(正增益准则)

4. 准则四(不相交准则)

1.5 基本LK算法的改进

5. 准则五

6. 准则六

7. 准则七

8. 准则八

9. 准则九

2. LKH算法

3. Python实现


0. 前言

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

Extreme Algorithmsicon-default.png?t=N7T8http://www.seas.gwu.edu/~simhaweb/champalg/tsp/tsp.html

LKH (Keld Helsgaun)icon-default.png?t=N7T8http://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的博客-CSDN博客_lkh算法

LKH求解TSP_yang_rookie的博客-CSDN博客_lkh算法

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

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

  • 4
    点赞
  • 33
    收藏
    觉得还不错? 一键收藏
  • 4
    评论
LKH(Lin-Kernighan)算法是一种用于求解TSP问题的优化算法,本质上是一种邻域搜索算法。下面是一个基于Python实现LKH算法的简单示例: ```python import random # 构造随机初始路径 def random_path(n): path = list(range(1, n+1)) random.shuffle(path) path.append(path[0]) return path # 计算路径长度 def path_length(path, dist): return sum(dist[path[i-1]-1][path[i]-1] for i in range(1, len(path))) # LK算法 def LK(dist, path, max_iters=10000): n = len(path)-1 for i in range(max_iters): improved = False for a in range(1, n+1): for b in range(a+1, n+1): # 1. 构造交换后的新路径 new_path = path[:a] + path[a:b+1][::-1] + path[b+1:] # 2. 计算新路径长度 new_length = path_length(new_path, dist) # 3. 判断是否接受新路径 if new_length < path_length(path, dist): path = new_path improved = True break if improved: break if not improved: break return path # 示例用法 if __name__ == '__main__': # 构造距离矩阵 dist = [[0, 3, 4, 2], [3, 0, 4, 3], [4, 4, 0, 2], [2, 3, 2, 0]] # 构造随机初始路径 path = random_path(len(dist)) # LK算法优化路径 path = LK(dist, path) # 输出结果 print(path) print(path_length(path, dist)) ``` 上述代码中,`random_path`函数用于构造随机初始路径,`path_length`函数用于计算路径长度。`LK`函数是LKH算法的核心实现,其中`max_iters`参数指定最大迭代次数,可以根据实际情况进行调整。最后,我们可以通过调用`LK`函数来使用LKH算法求解TSP问题。
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值