はじめに

本投稿は，初学者がインプットしたものをどうにかアウトプットの形にしたいという思いで投稿しています．そのため，やや間違っている部分もあるかと思います．
ぜひ，間違っていましたら，教えてください！

1. そもそも巡回セールスマン問題(TSP)って？

巡回セールスマン問題とは，複数の決められた地点をセールスマンが移動する際，どの順番が最短距離になるのか？という問題です．この問題，地点が増えれば増えるほど，計算量が増加します．(例えば，都市数をnとした時，考えられる解はn!/2nになります！）
そのため，解くのに非常に時間がかかるため，多項式時間に解けない「NP問題」と同等並に難しい「NP困難」に分類されています．

2. どんな解法があるの？

この問題を解くにあたって，大きく3つの解法があります．

厳密数値解法
ヒューリスティック解法
深層学習，強化学習を用いた方法

本投稿では，この3つの中のヒューリスティック解法の一種である，遺伝的アルゴリズムで解こうと思います．

3. 遺伝的アルゴリズムとは？

遺伝的アルゴリズムとは，生物の進化の仕組みを模した最適化アルゴリズムのことです．優れた親同士の特徴を子に受け継がせることで，より良い個体を作るというイメージです．(「テラフォーマーズ　ニュートン一族」で調べてもらえば，より分かりやすいかもしれません(わかる人にはわかる...))

流れとしては，

1. 初期個体の生成
2. 適応度の評価
3. 選択
4. 交叉
5. 突然変異
6. 世代交代

これの2~6を,「収束する or ある世代数」まで繰り返します．

3.1. 初期個体の生成

第1世代の個体を生成します．この時に作られた個体数が，以降の個体数になります．そのため個体数が多いと，計算時間は長くなりますが，その分，精度が高くなる傾向にあります．

3.2. 適応度の評価

生成された個体を評価します．今回の問題は，「移動距離を短くする」が目標なので，移動距離を評価します．

3.3. 選択

3.2での結果を用いて，交叉する親個体を選択します．いくつかの手法があり，

ルーレット選択
エリート選択
トーナメント選択

などがあります．(どの方法にどんな特徴があるかは結果を見ながら今後比較したいところ！)
今回は「ルーレット選択」を採用します．
※ルーレット選択とは，3.2の結果から確率的に交叉する親個体を選ぶ方法です．

3.4. 交叉

親個体の特徴を受け継ぐ子個体を生成します．この方法もいくつかあり，

部分写像交叉(PMX)
順序交叉(OX)
循環交叉(CX)
サブツア交換交叉
辺再組み合わせ交叉(ER)

等が存在します．(どの方法にどんな特徴があるかは結果を見ながら今後比較したいところ！)
今回は「順序交叉(OX)」を採用します．
※順序交叉(OX)とは，片方の親からは一部をそのまま受け継ぎ，残りの部分は相対的な順を保存しつつ受け継ぐ方法です．

3.5. 突然変異

親の個体を無視して，新しく個体を生成します．これは，親の特徴を受け継ぐ子のみを生成し続けると，局所解に入る恐れがあるからです．

3.6. 世代交代

以上の結果で生成した個体を次の親世代とし，3.2~3.5を繰り返します．

4. 実装結果

では，実際にTSPを解いていきます．今回は，「TSPLIB」というサイトから問題を選びました．このサイトでは，WPの初期配置とその時の最短経路が記載されているため，比較しやすいです．
問題設定としては，下の図のようになっています．

・「berlin52」

そして，最適解は以下のようになっています．

パラメータは以下のようにしています．

世代数：100000
親世代をそのまま子世代に受け継がせる数：10
突然変異する確率：1%
個体数：25
試行回数：5

最短経路の距離が7544.37であるので，(TSPLIBの評価関数の値とは異なるので，要検討ではあります)これとの比較を行います．

こんな感じで，少しずつ最短距離に近づいている様子が確認できます．

以下，5回試行した場合の結果を示します．

＜1回目＞

＜２回目＞

＜３回目＞

＜４回目＞

＜５回目＞

これらの試行が世代数によってどのように変化したか，そして，最短経路とどれ程の差異があるのかを比較したグラフが以下の図です．

これを見ると，10000世代ぐらいまでは大きく変化をしていますが，それ以降は少しずつしか変化していないことが見て取れます．
また，最短経路にはいずれも到達していないことはわかります．
このため，選択方法や交叉方法を別のものにすることを検討する必要があると考えています．また，突然変異のフェーズにおいて，その世代の最適解が失われるようなプログラムになっていることも原因の一つだと考えています．
あと，収束条件も入れる必要があるは思っています．(全体で70分程度かかっているので...)

5. 今後の課題

今後の課題としては

TSPLIBの評価関数と同じ評価方法を用いる
選択方法，交叉方法を変更する
世代数などのパラメータの最適解を見つける
ヒューリスティック解法だけではなく，今話題のAIを用いた方法も試してみる

等を行いたいです．

6. 使用プログラム

以下，今回用いたコードです．

#GA_TSP
import numpy as np
import random

#遺伝子を生成
def generate_init_genes(num_individual, num_cities):
    genes = np.zeros((num_individual, num_cities), dtype=np.int16)
    for i in range(num_individual):
        genes[i,] = random.sample(range(num_cities), k=num_cities)
    return genes

#都市間の距離についてのテーブル
def generate_dist_table(positions):
    path_table = np.zeros((len(positions), len(positions)))
    for i in range(len(positions)):
        for j in range(len(positions)):
            path_table[i,j] = np.linalg.norm(positions.loc[i] - positions.loc[j])
    return path_table

#テーブルを参照して、遺伝子に対応させた経路の和を求める
def sum_path2(path_table, gene):
    sum = 0.
    for i in range(len(gene)-1):
        sum += path_table[int(gene[i]),int(gene[i+1])]
    return sum

#経路の和について、一世代計算
def genes_path2(path_table, genes):
    pathlength_vec = np.zeros(len(genes))
    for i in range(len(genes)):
        pathlength_vec[i] = sum_path2(path_table, genes[i])
    return pathlength_vec

#テーブルを参照して、遺伝子に対応させた経路の和を求める
def sum_path2_lap(path_table, gene):
    sum = 0.
    for i in range(len(gene)-1):
        sum += path_table[int(gene[i]),int(gene[i+1])]
    sum += path_table[gene[-1], gene[0]]
    return sum

#経路の和について、一世代計算
def genes_path2_lap(path_table, genes):
    pathlength_vec = np.zeros(len(genes))
    for i in range(len(genes)):
        pathlength_vec[i] = sum_path2_lap(path_table, genes[i])
    return pathlength_vec

#ルーレット
def generate_roulette(fitness_vec):
    total = np.sum(fitness_vec)
    roulette = np.zeros(len(fitness_vec))
    for i in range(len(fitness_vec)):
        roulette[i] = fitness_vec[i]/total
    return roulette

#ルーレットをもとに交叉をする個体を選ぶ
def roulette_choice(fitness_vec):
    roulette = generate_roulette(fitness_vec)
    choiced = np.random.choice(len(roulette), 2, replace=False, p=roulette)
    return choiced

def order_crossover(parent1, parent2):
    num = len(parent1)

    # ランダムな交叉範囲を決定
    cross_range = random.sample(range(num), 2)
    start, end = min(cross_range), max(cross_range)

    # 子供の初期化
    child1 = [-1] * num
    child2 = [-1] * num

    # 交叉範囲の遺伝子をコピー
    child1[start:end] = parent1[start:end]
    child2[start:end] = parent2[start:end]

    # 親の残りの遺伝子を使って子供の残り部分を埋める
    def fill_remaining(child, parent):
        current_pos = end  # 交叉範囲の終了位置からスタート
        for gene in parent:
            if gene not in child:  # 重複しないように配置
                if current_pos >= num:
                    current_pos = 0
                child[current_pos] = gene
                current_pos += 1
    
    fill_remaining(child1, parent2)
    fill_remaining(child2, parent1)

    return child1, child2

def translocation_mutation2(genes, p_value):
    mutated_genes = genes
    for i in range(len(genes)):
        mutation_flg = np.random.choice(2, 1, p = [1-p_value, p_value])
        if mutation_flg == 1:
            mutation_value = np.random.choice(genes[i], 2, replace = False)
            mutation_position1 = np.where(genes[i] == mutation_value[0])
            mutation_position2 = np.where(genes[i] == mutation_value[1])
            mutated_genes[i][mutation_position1] = mutation_value[1]
            mutated_genes[i][mutation_position2] = mutation_value[0]
    return mutated_genes

def translocation_mutation_multi(genes, p_value, num_points):
    mutated_genes = genes
    for i in range(len(genes)):
        mutation_flg = np.random.choice(2, 1, p = [1-p_value, p_value])
        if mutation_flg == 1:
            print(i)
            # ランダムにnum_points個の位置を選択し、その部分をシャッフル
            mutation_indices = np.random.choice(len(genes[i]), num_points, replace=False)
            #print(mutation_indices)
            mutation_values = mutated_genes[i][mutation_indices]
            np.random.shuffle(mutation_values)
            # 元の遺伝子配列にシャッフル後の値を代入
            mutated_genes[i][mutation_indices] = mutation_values

    return mutated_genes

import matplotlib.pyplot as plt
import pandas as pd
import GA_TSP
import numpy as np
import time
import os

def show_route(cities, gene, distance, title="Optimal Route"):
    plt.scatter(cities[:, 0], cities[:, 1], c='blue', marker='o')
    
    # 開始地点に戻るように、経路の最後に開始地点（gene[0]）を追加
    gene = np.append(gene, gene[0])
    
    for i in range(len(gene) - 1):
        if i == 0:
            plt.text(cities[int(gene[i])][0], cities[int(gene[i])][1], "start", fontsize=10, color="green")
        else:
            plt.text(cities[int(gene[i])][0], cities[int(gene[i])][1], str(i), fontsize=8, color="blue")
        plt.plot([cities[int(gene[i])][0], cities[int(gene[i+1])][0]], 
                 [cities[int(gene[i])][1], cities[int(gene[i+1])][1]], 'k-')
    plt.text(cities[int(gene[-1])][0], cities[int(gene[-1])][1], "goal", fontsize=10, color="red")   
    plt.title(title)
    plt.xlabel("X")
    plt.ylabel("Y")

# 計算開始
start_time = time.time()

fname = 'berlin52.tsp/berlin52.csv'
wp_data = pd.read_csv(fname, header=None, index_col=0, encoding="cp932")
wp_data = wp_data.reset_index(drop=True)

# 初期配置のプロット
plt.figure()
plt.scatter(wp_data[1], wp_data[2], c='blue', marker='o', label="Waypoints")
plt.xlabel("X")
plt.ylabel("Y")
plt.legend()
plt.show()

# パラメータ設定
generation = 100000
elite = 10
p_mutation = 0.01
init_individual = 25
trials = 5  # 試行回数

# 各試行のベスト距離の推移と最適経路を格納するリスト
all_trials_best_distances = []
all_trials_best_genes = []  # 各試行の最適経路の遺伝子配列
all_trials_best_distance_values = []  # 各試行の最小距離

# 結果保存のためのフォルダーを作成
output_dir = "GA_TSP_Results"
os.makedirs(output_dir, exist_ok=True)

for trial in range(trials):
    max_fit = 0
    trial_folder = os.path.join(output_dir, f"Trial_{trial+1}")
    os.makedirs(trial_folder, exist_ok=True)  # 試行ごとのフォルダー作成

    # 初期個体の生成
    genes = GA_TSP.generate_init_genes(init_individual, wp_data.shape[0])
    dist_table = GA_TSP.generate_dist_table(wp_data)
    paths = GA_TSP.genes_path2(dist_table, genes)
    inverse_path = 1 / paths

    # ルーレット選択の初期化
    GA_TSP.generate_roulette(inverse_path)

    # 実行部
    fitness_vec = 1 / (GA_TSP.genes_path2(dist_table, genes))
    best_distances = []

    for i in range(generation):
        print(i)
        children = np.zeros(np.shape(genes))
        for j in range(int((init_individual - elite) / 2 + 1)):
            parents_indices = GA_TSP.roulette_choice(fitness_vec)
            children[2 * j], children[2 * j + 1] = GA_TSP.order_crossover(genes[parents_indices[0]], genes[parents_indices[1]])

        # エリート選択
        for j in range(init_individual - elite, init_individual):
            children[j] = genes[np.argsort(GA_TSP.genes_path2(dist_table, genes))[j - init_individual + elite]]

        genes = children
        fitness_vec = 1 / (GA_TSP.genes_path2(dist_table, genes))
        
        # 開始地点に戻る経路の距離を計算
        best_gene = genes[np.argmax(fitness_vec)]
        best_distance = GA_TSP.sum_path2(dist_table, best_gene) + dist_table[int(best_gene[-1]), int(best_gene[0])]
        best_distances.append(best_distance)
        
        # 過去最高の適応度を記録した場合、経路を表示し保存
        if max(fitness_vec) > max_fit:
            max_fit = max(fitness_vec)
            plt.clf()  # グラフをクリア
            plt.scatter(wp_data[1], wp_data[2], c='blue', marker='o')
            show_route(wp_data.values, genes[np.argmax(fitness_vec)], best_distance)

            # テキスト表示
            mid_x, mid_y = np.mean(wp_data[1]), np.mean(wp_data[2])  # グラフの中心を計算
            plt.text(mid_x, mid_y, f"Generation: {i}\nDistance: {best_distance:.2f}",
                     fontsize=12, color="black", ha='center', bbox=dict(facecolor='white', alpha=0.5))
            
            save_path = os.path.join(trial_folder, f"Generation_{i}_Route.png")
            plt.savefig(save_path)  # 画像を保存
            plt.pause(1)  # 1秒の待機

    # 各試行のベスト距離の推移と最短経路を保存
    all_trials_best_distances.append(best_distances)
    all_trials_best_genes.append(best_gene)
    all_trials_best_distance_values.append(best_distance)

    print(f"Trial {trial+1} completed")

# 最短距離の推移をプロット
plt.figure()
for i, best_distances in enumerate(all_trials_best_distances):
    plt.plot(range(generation), best_distances, label=f"Trial {i+1}")

# y=7544.37のラインを追加
target_distance = 7544.37
plt.axhline(y=target_distance, color='red', linestyle='--', label=f'y = {target_distance}')

# 計算時間の修了
end_time = time.time()
calcuration_time = end_time - start_time
print(f"{calcuration_time:.4f} seconds")

plt.xlabel("Generation")
plt.ylabel("Best Distance")
plt.title("Best Distance vs. Generation Across Multiple Trials in GA")
plt.legend()
save_path = os.path.join(output_dir, f"Trial_Optimal_Route_vs_Optimal_Route.png")
plt.savefig(save_path)

# 各試行の最短経路を個別のウィンドウに表示し保存
for trial, (gene, distance) in enumerate(zip(all_trials_best_genes, all_trials_best_distance_values)):
    plt.figure()  # 新しいウィンドウを作成
    show_route(wp_data.values, gene, distance, title=f"Optimal Route for Trial {trial+1}")
    # テキスト表示
    mid_x, mid_y = np.mean(wp_data[1]), np.mean(wp_data[2])  # グラフの中心を計算
    plt.text(mid_x, mid_y, f"Trial: {trial + 1}\nDistance: {distance:.2f}",
                fontsize=12, color="black", ha='center', bbox=dict(facecolor='white', alpha=0.5))
    save_path = os.path.join(output_dir, f"Trial_{trial + 1}_Optimal_Route.png")
    plt.savefig(save_path)  # 試行ごとの最適経路画像を保存
plt.show()

7. 参考文献

以下，TSP，GAを勉強する上で参考にした資料です．

「pythonで遺伝的アルゴリズム（GA）を実装して巡回セールスマン問題（TSP）をとく」をやってみた
https://qiita.com/spwimdar/items/e46487e25f35953c7cfe
遺伝的アルゴリズムの基礎　大阪市立大学
https://www.econ.osaka-cu.ac.jp/~hashimo/Genetic.pdf
TSPLIB Heidelberg University
http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/
交叉と突然変異
https://www.obitko.com/tutorials/genetic-algorithms/japanese/crossover-mutation.php