使用模拟退火和遗传算法求解 TSP 问题

import time
import io
import numpy
import random
from matplotlib import pyplot

def get_coords(s):
    with io.open(s, 'r') as f:
        coords = []
        coords_flag = 0
        for s in f.readlines():
            if s.strip() == 'EOF':
                coords_flag = 0
            if coords_flag:
                t = s.split()
                coords.append([int(t[1]), int(t[2])])
            if s.strip() == 'NODE_COORD_SECTION':
                coords_flag = 1
        return coords


def cal_length(solution, adjmap):
    length = 0
    for i in range(len(solution)):
        length += int(adjmap[solution[i-1]][solution[i]])
    return length


def tsp(method, test_time=10, coords=get_coords('eil101.tsp')):
    adjmap = [[numpy.linalg.norm(numpy.array(ci)-numpy.array(cj))
               for cj in coords]for ci in coords]
    length = []
    solution = []
    run_time = []
    ans = []
    for i in range(test_time):
        start = time.perf_counter()
        solution1, length1 = method(adjmap)
        run_time.append(time.perf_counter()-start)
        ans.append(length1[-1])
        if i == 0 or length[-1] > length1[-1]:
            length = length1
            solution = solution1

    print('以下是各次运行时间（s）及求解时间的均值、方差、标准差：')
    print(run_time)
    print(numpy.mean(run_time))
    print(numpy.var(run_time))
    print(numpy.std(run_time))

    print('以下是各次运行求得解的大小及其均值、方差、标准差：')
    print(ans)
    print(numpy.mean(ans))
    print(numpy.var(ans))
    print(numpy.std(ans))
    print(solution)

    print('以下是所求最优解及其可视化：')
    print(length[-1])
    x = [coords[i][0]for i in solution]
    y = [coords[i][1]for i in solution]
    x.append(x[0])
    y.append(y[0])
    pyplot.plot(x, y)
    pyplot.show()

    print('以下是最优解的收敛过程：')
    pyplot.plot(length)
    pyplot.show()

def variation(solution):
    def swap_pos(solution):
        new_s = solution.copy()
        pos = random.sample(range(len(new_s)), 2)
        pos.sort()
        new_s[pos[0]], new_s[pos[1]] = new_s[pos[1]], new_s[pos[0]]
        return new_s

    def reverse_pos(solution):
        new_s = solution.copy()
        pos = random.sample(range(len(new_s)), 2)
        pos.sort()
        new_s1 = new_s[0:pos[0]:1]
        new_s2 = new_s[pos[0]:pos[1]:1]
        new_s3 = new_s[pos[1]:len(new_s):1]
        return new_s1+list(reversed(new_s2))+new_s3

    def flip_pos(solution):
        new_s = solution.copy()
        pos = random.sample(range(len(new_s)), 2)
        pos.sort()
        new_s = new_s[pos[1]:len(new_s):1] + \
            new_s[pos[0]:pos[1]:1]+new_s[0:pos[0]:1]
        return new_s

    def rotate_pos(solution):
        new_s = solution.copy()
        pos = random.sample(range(len(new_s)), 2)
        pos.sort()
        new_s = new_s[pos[0]:pos[1]:1] + \
            new_s[0:pos[0]:1]+new_s[pos[1]:len(new_s):1]
        return new_s
    shuffer = random.choice([swap_pos, reverse_pos, flip_pos, rotate_pos])
    return shuffer(solution)


def local_search(adjmap, num_epoch=100000):
    # 构造一组初始解
    solution = [i for i in range(len(adjmap))]
    length = [cal_length(solution, adjmap)]
    for i in range(1, num_epoch+1):
        new_s = variation(solution)
        new_l = cal_length(new_s, adjmap)
        if new_l < length[-1]:
            solution = new_s
            length.append(new_l)
        else:
            length.append(length[-1])
    return solution, length


tsp(local_search)

def simulated_annealing(adjmap, t_beg=100, t_end=0.01, q=0.99, l=100):
    # 构造一组初始解
    solution = [i for i in range(len(adjmap))]
    length = [cal_length(solution, adjmap)]
    t = t_beg
    while t > t_end:
        for i in range(l):
            new_s = variation(solution)
            new_l = cal_length(new_s, adjmap)
            pre_l = cal_length(solution, adjmap)
            df = new_l-pre_l
            if df < 0 or numpy.exp(-df / t) > random.random():  # 退火，一定概率下接受差解
                solution = new_s
        length.append(cal_length(solution, adjmap))
        t *= q
    return solution, length


tsp(simulated_annealing)

def cross(p1, p2):
    a = p1.copy()
    b = p2.copy()

    pos = random.sample(range(len(p1)), 2)
    pos.sort()

    # print pos[0],end
    # 建立映射关系
    cross_map = {}
    is_exist = False
    # 初步映射
    for i in range(pos[0], pos[1]):
        if a[i] not in cross_map.keys():
            cross_map[a[i]] = []
        if b[i] not in cross_map.keys():
            cross_map[b[i]] = []

        cross_map[a[i]].append(b[i])
        cross_map[b[i]].append(a[i])

    # 处理传递映射 如1:[6],6:[3,1]-->1:[6,3,1],6:[3,1]
    # 计算子串中元素出现的个数，个数为2，则该元素为传递的中间结点，如如1:[6],6:[3,1],‘6’出现的次数为2
    appear_times = {}
    for i in range(pos[0], pos[1]):
        if a[i] not in appear_times.keys():
            appear_times[a[i]] = 0
        if b[i] not in appear_times.keys():
            appear_times[b[i]] = 0

        appear_times[a[i]] += 1
        appear_times[b[i]] += 1

        if a[i] == b[i]:
            appear_times[a[i]] -= 1

    for k, v in appear_times.items():
        if v == 2:
            values = cross_map[k]
            for key in values:
                cross_map[key].extend(values)
                cross_map[key].append(k)
                cross_map[key].remove(key)
                cross_map[key] = list(set(cross_map[key]))

    # 使用映射关系交叉
    # 先映射选中的子串
    temp = a[pos[0]:pos[1]].copy()
    a[pos[0]:pos[1]] = b[pos[0]:pos[1]]
    b[pos[0]:pos[1]] = temp

    # 根据映射规则映射剩下的子串
    seg_a = a[pos[0]:pos[1]]
    seg_b = b[pos[0]:pos[1]]

    remain = list(range(pos[0]))
    remain.extend(range(pos[1], len(a)))

    for i in remain:
        keys = cross_map.keys()
        if a[i] in keys:
            for fi in cross_map[a[i]]:
                if fi not in seg_a:
                    a[i] = fi
                    break
        if b[i] in keys:
            for fi in cross_map[b[i]]:
                if fi not in seg_b:
                    b[i] = fi
                    break
    return [a, b]


def genetic_algorithm(adjmap, num_epoch=25000, chose_num=4, up_bound=10, p_c=0.9, p_m=0.4):
    # 构造初始种群
    w = [i for i in range(len(adjmap))]
    population = []
    for i in range(up_bound):
        population.append(w.copy())
        w = variation(w)
    length = []
    for i in range(num_epoch):
        # 适应者繁殖
        w = [cal_length(p, adjmap) for p in population]
        sum = 0
        for ww in w:
            sum += ww
        chose = random.choices(population, weights=[
                               numpy.exp(-ww/sum) for ww in w], k=chose_num)
        for j in range(1, chose_num, 2):
            parent = [chose[j-1], chose[j]]
            if random.random() < p_c:
                parent = cross(parent[0], parent[1])
            for ch in parent:
                if random.random() < p_m:
                    ch = variation(ch)
                population.append(ch)
        # 不适者淘汰
        population.sort(key=lambda sol: cal_length(sol, adjmap))
        while len(population) > up_bound:
            population.pop()
        length.append(cal_length(population[0], adjmap))
    return population[0], length


tsp(genetic_algorithm)