pythonで遺伝的アルゴリズム(GA)を実装して巡回セールスマン問題(TSP)をとく

遺伝的アルゴリズムについての説明と実装した結果について紹介します。
具体例として解く問題を巡回セールスマン問題にしています。

目次

  1. KUB
  2. 巡回セールスマン問題(TSP)について
  3. 遺伝的アルゴリズム(GA)について
  4. 遺伝的アルゴリズムの実装
  5. まとめ

1. KUB

皆さんはKUBという略語をご存知でしょうか。
KUBは私が今作った言葉で、
こんな胡散臭いブログ参考になるか!/参考にしたくない!」の略です。
そんな方のためにもっと正しい情報にアクセスできるような情報源を
紹介しようと思います。

TSP Tutorial in Python by Peter Norvig
jupyter形式のTSPのわかりやすい解説。
遺伝的アルゴリズムこそありませんが、
総当り法やNN法などを利用しています。
言語は英語ですがコードや図を見ればなんとなくわかると思います。

DEAP (Github) (documentation)
遺伝的アルゴリズムを簡単に実装できる進化計算用のフレームワーク。
解説(英語)が公式できちんと出ている上に、日本語のやってみた記事が
いくつもヒットするので敷居は低めだと思います。

以降の記事でもKUBのコーナーを設けるつもりなのでよしなに。

2. 巡回セールスマン問題(TSP)

以前の記事でTSP(Time Stretched Pulse)信号を紹介しましたが、
これとはまったく関係がありません。

巡回セールスマン問題とは、組み合わせ最適化問題のひとつです。
組み合わせ最適化問題を平たく言うと、とある条件にあう組み合わせの内で、
設定した得点(コスト)が最大(最小)になる組み合わせを探す、
という問題です。
他の組み合わせ最適化問題だとナップサック問題が有名です。

具体的には、都市が複数与えられて、その都市をすべて一回ずつ通る、
一番経路が短い巡り方はどれか?という問題です(Fig. 1)。

Fig. 1 TSPの簡単な説明

このTSPが一体どうして幅広く問題として用いられるかというと、
シンプルでわかりやすい割に計算量が莫大に増えるからです。

Fig. 2 都市の数の階乗で増える距離の組み合わせ

Fig. 2 に都市の数と距離の組み合わせの数のイメージを示します。
都市が2つの場合は2通りの答えしか出ません。
都市が3つだと6通り、4つだと24通り、5つだと120通り、
というように都市の数をnとおくと距離の組み合わせはn!通りです。

まだ120くらいならいいですが、20都市でやろうとすると、
2432902000000000000通り距離の組み合わせが存在します(多分)。
このように数え上げると莫大な問題は逐一計算すると、
解くのに非現実的な時間がかかります。
このような問題を解くために様々なアルゴリズムが考案されています。
(参考:「『フカシギの数え方』 おねえさんといっしょ! みんなで数えてみよう!」
日本科学未来館の公式動画YouTubeページ
ぶっ飛んだ角度からアルゴリズムの大切さがわかります。)

3. 遺伝的アルゴリズム(GA)

依頼主はいつものGA

遺伝的アルゴリズム(Genetic Algorithm)とは、
生物の進化のメカニズムを取り入れた最適化アルゴリズムのひとつです。

最適化というのは決められたルールによって計算される得点を最大にする
(もしくはコストを最小にする)ことを指します。

アルゴリズム特定のタスクを達成するための手順です。

生物の進化のメカニズムはダーウィンの進化論にあたり、
その環境に適応できなかった個体が淘汰され(自然淘汰)、
残った個体間でそれらの遺伝子の交叉が行われて、
次世代に引き継がれるというものです。

生物は染色体という生体物質の中に、その個体の設計図ともなる遺伝子情報
をもっています。詳しくはMSDマニュアル家庭版「遺伝子と染色体」へ。

例えば、生物が[暑くても平気][寒くても平気][短期間なら両方平気]の
三種類のいずれかをもっているとします。
この生物たちは現状ではバランスよく存在していますが、
あるとき気候変動が起こり、暑い環境になったとします。
そうすると暑さに耐えられる個体が生き残るようになります。(Fig. 3)

Fig. 3 自然淘汰のイメージ

これが自然淘汰です。
この場合、対応すべき環境は暑さであり、暑さ耐性を持つことで生き残れます。
生き残った個体間で交配が行われ、暑さ耐性を持つ個体が増えます。(Fig. 4)

Fig. 4 環境に適応して生き残った個体の特徴を引き継いだ次世代が誕生

この流れ(環境に適応できているかの判断→淘汰→交叉)を繰り返すことで、
この生物はみんな環境に適応することができます。
しかし、この戦略には弱点があります。仮にしばらく暑い期間が続いたのちに、
急激に寒くなった場合を考えてみます。考えるまでもなく全滅です。(Fig. 5)

Fig. 5 世代が進むと集団が同一になり環境変化に脆弱になるイメージ

生物は突然変異という機能によってこの弱点を克服しています。
突然変異とは平たく言うと「たまに違うやつが現れる」ことです。(Fig. 6)

Fig. 6 突然変異で全滅を免れる生物たち

遺伝的アルゴリズムは以上にあげたような、
環境への適応度を求める評価、生き残る個体を決める選択
生き残った個体で遺伝情報を交換する交叉
そして多様性を担保するための突然変異
これら4つのプロセスを繰り返すことで構成されています。

4. 遺伝的アルゴリズムの実装

遺伝的アルゴリズムのフローチャートはFig. 7のようになっています。
今回は解く問題を巡回セールスマン問題に限定しているので、
より具体的な内容も併記しました。

Fig. 7 遺伝的アルゴリズムの一般的なフローチャート(左)
と今回のTSPのための補足(右)

TSPの設定

まずTSPの問題設定のために指定した数だけ都市の位置(x, y)を
ランダムで生成する関数を実装します。

今回はx, yともに0~1の範囲内に都市があるとしています。

入力は都市数、出力は都市の位置を指す都市数×2の行列です。

実装例: generate_rand_cities

def generate_rand_cities(num_cities):
    positions = np.zeros((points, 2))
    for i in range(points):
        positions[i, 0] = random.random()
        positions[i, 1] = random.random()
    return positions

使用例: generate_rand_cities

generate_rand_cities(10)

都市の数を10とした場合です。

使用例: generate_rand_cities  実行結果

array([[0.07730533, 0.27874579],
       [0.58098837, 0.54837502],
       [0.59160204, 0.47646458],
       [0.51356026, 0.85514733],
       [0.4312847 , 0.39608552],
       [0.13582582, 0.96681264],
       [0.84180448, 0.38115363],
       [0.93112806, 0.31603043],
       [0.96833881, 0.38714716],
       [0.14744214, 0.32422133]])

遺伝子の構成と初期化

Fig. 8 設定した遺伝子の構成

Fig. 8に示したような遺伝子の構成を用います。
一個体の遺伝子は、巡る都市が順々に並んでいるベクトルです。
したがってその長さは都市の数と等しくなります。
例えば、都市が5つ存在している場合は、
[0, 1, 2, 3, 4]
[2, 4, 1, 0, 3]
などのような個体が考えられます。

また一世代の内に複数個体を作るのでその個体数をm、
都市の数をnとすると一世代の遺伝子はm×nの行列で表せます。

都市の数と個体数を入力して初期個体をランダムで生成して
行列の形で出力する関数を実装します。

実装例: generate_init_genes

def generate_init_genes(num_indivisual, num_cities):
    genes = np.zeros((num_indivisual, num_cities), dtype=np.int16)
    for i in range(num_indivisual):
        genes[i,] = random.sample(range(num_cities), k=num_cities)
    return genes

使用例: generate_init_genes

generate_init_genes(15, 10)

一世代の個体数を15、都市の数を10にした例。

使用例: generate_init_genes 実行結果

array([[5, 9, 8, 0, 4, 1, 3, 2, 7, 6],
       [6, 1, 8, 9, 0, 4, 3, 5, 2, 7],
       [5, 9, 7, 4, 3, 6, 8, 2, 0, 1],
       [7, 8, 1, 3, 0, 2, 9, 4, 5, 6],
       [5, 6, 3, 0, 8, 1, 2, 7, 4, 9],
       [5, 2, 0, 7, 9, 1, 4, 6, 8, 3],
       [2, 4, 1, 3, 5, 7, 6, 0, 9, 8],
       [7, 3, 9, 0, 4, 1, 8, 5, 6, 2],
       [0, 9, 5, 2, 6, 8, 4, 1, 7, 3],
       [4, 5, 0, 9, 7, 2, 1, 3, 8, 6],
       [5, 7, 3, 4, 8, 1, 0, 2, 9, 6],
       [5, 9, 4, 7, 3, 2, 8, 0, 6, 1],
       [5, 8, 4, 6, 7, 1, 9, 2, 0, 3],
       [0, 1, 6, 8, 3, 9, 4, 7, 2, 5],
       [2, 4, 7, 6, 9, 5, 1, 8, 0, 3]], dtype=int16)

各要素に格納されている数字は巡る都市を指しています。
したがってint型であることを明示的にしています。
またここで登場している数字はgenerate_rand_citiesの行数と
対応させて使用します。

今回はTSPに即した遺伝子の構成であり、
他の構成の仕方としては複数同じものが選ばれるもの、
ナップザック問題で使われる対立遺伝子が0と1で表現されるもの
などさまざまあります。

評価関数の設定

ここまでで都市と、その巡る順番を指定することが
できるようになりました。

次はその遺伝子がどれほど環境に適合しているか?
今回のTSPに即して言えばどれほど経路が短いか?を評価します。
といっても距離を測るだけなのであまり難しくはないです。

賢い方法としては2点間の距離のテーブルを作って、
順序に合わせて参照する方法があります。
が、今回は頭良くない方法でいきます。
具体的にはその都度、都市を参照して距離を計算します。

まずは一個体の距離の総和を求める関数から実装します。
入力は都市の位置と遺伝子、出力は距離の値です。

実装例: sum_path

def sum_path(positions, gene):
    sum = 0.
    for i in range(len(positions)-1):
        sum += np.linalg.norm(positions[gene[i]]- 
                              positions[gene[i+1]])
    return sum

2点間の距離を求めるのにnumpy.linalg.normを用いました。

使用例: sum_path

cities = generate_rand_cities(10)
genes = generate_init_genes(15, 10)
for i in range(15):
    print(sum_path(cities, genes[i]))

同じく個体数15、都市数10。

使用例: sum_path 実行結果

3.255120729540855
4.101759745461325
3.707831808175001
5.292216946473589
4.499663586254492
5.175813297976937
4.801956481827147
4.374668651228482
4.466671685532946
3.8037980851237516
4.4033551408241
4.818583924144366
4.731557533385495
4.648770612139303
4.160335620768928

使用例では個体数の分だけfor文で計算して個々の遺伝子に対して
経路の長さを求めています。
for文書くのも野暮ったいしベクトルとしてまとめて出力したいので
さらに一世代分まるごと経路の長さを求めて出力する関数を実装します。

実装例: genes_path

def genes_path(genes, cities):
    pathlength_vec = np.zeros(len(genes))
    for i in range(len(genes)):
        indices = genes[i]
        pathlength_vec[i] = sum_path(cities, indices)
    return pathlength_vec

使用例: genes_path

cities = generate_rand_cities(10)
genes = generate_init_genes(15, 10)

genes_path(genes, cities)

同じく個体数15、都市数10。

使用例: genes_path 実行結果

array([5.08913991, 5.3178538 , 4.78949816, 4.64098699, 6.03794245,
       5.74928871, 3.9795153 , 4.66686211, 5.06521305, 5.72204159,
       4.64535672, 5.09412403, 5.03756046, 4.89669194, 4.89986321])

選択(ルーレット、エリート)

選択の内でも交叉する両親を選ぶものと、
ダイレクトに子孫として残す個体を選ぶものの2種類存在します。
前者にはルーレット選択やランキング選択などがあり、
後者にはエリート選択があります(Fig. 9)。

Fig. 9 選択の違いのイメージ

ルーレット選択もランキング選択も適応度が高い個体が選ばれる
確率が高なるような選択方法です。

Fig. 10 ルーレット選択の具体例とイメージ

ルーレット選択は適応度に応じて選択される確率が変わります。

実装例: generate_roulette

def genarate_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

使用例: generate_roulette

fitness = np.array([20,50,30])
generate_roulette(fitness)

Fig. 10 の例と同じように適応度が20、50,30の3つの個体があったとき、
それぞれが選択される確率を出力します。

使用例: generate_roulette 実行結果

array([0.2, 0.5, 0.3])

今回は移動距離をもとに適応度を計算しますが、
移動距離は短ければ短いほどいいのに対して、
適応度が高ければ高いほど選択されるため、少し細工が必要です。
その細工として逆数を利用した例とその結果を以下に示します。
(個体数: 3、都市数: 5 として)

cities = generate_rand_cities(5)
genes = generate_init_genes(3, 5)

path = genes_path(genes, cities)
inverse_path = 1/path
print("path length: "+str(path))
print("roulette table: "+str(generate_roulette(inverse_path)))

実行結果

path length: [2.76269891 2.89223783 1.45059135]
roulette table: [0.25908451 0.24748051 0.49343497]

道のり(path length)の小さい要素は、選択される確率(roulette table)が
高くなっていることを確認できます。

ルーレット、つまり選ばれる確率のテーブルができたので、
これを元に交叉を行う個体を選択する関数を実装します。

実装例: roulette_choice

def roulette_choice(fitness_vec):
    roulette = generate_roulette(fitness_vec)
    choiced = np.random.choice(len(roulette), 2, replace=True, p=roulette)
    return choiced

2行目でgenerate_rouletteを使ってルーレットを生成し、
またnumpy.random.choiceを使ってルーレットを実行しています。
第1引数が選ぶ数字のレンジで、rouletteの長さ、つまり個体数があてはまり、
第2引数は選ぶ数字の数、replace=は重複を許すか否かのオプションで、
p=でそれぞれが選ばれる確率を指定しています。

使用例: roulette_choice

cities = generate_rand_cities(10)
genes = generate_init_genes(6, 10)
fitness_vec = 1 / genes_path(genes, cities)
for i in range(3):
    print(roulette_choice(fitness_vec))

使用例: roulette_choice 実行結果

[2 0]
[1 1]
[4 3]

エリート選択は成績の良い個体をそのまま保存することを目的としているため
以下に示す交叉や突然変異の対象にはしません。

交叉(部分的交叉)

交叉にも遺伝子の形によって適した方法があったり、
あるいは同じ遺伝子の形に対しても様々な方法があったりします。
色々変えて適切なのはどれか、選べる環境がベストだと思います。

さて今回の遺伝子の形は、
都市のインデックスを巡る順番に格納しているというものです。
順番の意味を失うことなく、遺伝子が重複しない手法をとる必要があります。

そのような交叉方法として、循環交叉や順序交叉などが挙げられますが、
部分的交叉という手法を実装します。

部分的交叉の具体的な例を用いた手順をFig. 11に示します。
(画像自体は4分割されています。)

Fig. 11 部分的交叉の具体的な手順

実装例: partial_crossover

def partial_crossover(parent1, parent2):
    num = len(parent1)
    cross_point = random.randrange(1, num-1)
    child1 = parent1
    child2 = parent2
    for i in range(num - cross_point):
        target_index = cross_point + i
        
        target_value1 = parent1[target_index]
        target_value2 = parent2[target_index]
        exchange_index1 = np.where(parent1 == target_value2)
        exchange_index2 = np.where(parent2 == target_value1)

        child1[target_index] = target_value2
        child2[target_index] = target_value1
        child1[exchange_index1] = target_value1
        child2[exchange_index2] = target_value2
    return child1, child2

numpy.whereは指定した値が対象とする行列の中のどの位置にあるか
インデックスで返してくれるものです。 (Fig. 11の③の動作)

使用例: partial_crossover

genes = generate_init_genes(2, 10)
print("parent1: "+str(genes[0]))
print("parent2: "+str(genes[1]))
child = partial_crossover(genes[0], genes[1])
print("child1:  "+str(child[0]))
print("child2:  "+str(child[1]))

使用例: partial_crossover 実行結果

parent1: [1 7 3 8 2 6 9 0 4 5]
parent2: [0 2 1 7 9 3 5 4 6 8]
child1:  [1 7 3 9 2 0 5 4 6 8]
child2:  [4 2 1 7 5 3 8 6 0 9]

突然変異(転座)

突然変異にも様々な手法が存在します。
転座はランダムで選んだ2つの要素を交換するものです。

Fig. 12 転座の具体的なイメージ

実装例: translocation_mutation

def translocation_mutation(genes, num_mutation, p_value):
    mutated_genes = genes
    for i in range(num_mutation):
        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

使用例: translocation_mutation

genes = generate_init_genes(5, 10)
print("before")
print(genes)
print("after")
print(translocation_mutation(genes, 3, 0.7))

個体数5、都市数10、突然変異がおこるかもしれない個体数3、
突然変異の確率を0.7に設定しています。

実装例: translocation_mutation 実行結果

before
[[5 6 8 2 4 0 1 3 7 9]
 [8 6 1 0 5 3 9 4 2 7]
 [3 2 8 9 7 4 0 1 6 5]
 [9 4 6 1 2 7 5 8 3 0]
 [9 7 8 0 1 4 5 3 6 2]]
after
[[5 6 8 0 4 2 1 3 7 9]
 [8 6 1 0 5 3 9 4 2 7]
 [3 4 8 9 7 2 0 1 6 5]
 [9 4 6 1 2 7 5 8 3 0]
 [9 7 8 0 1 4 5 3 6 2]]

上から3番目までは変化があったりなかったり。
繰り返し実行してみるとわかりやすいかと思います。

可視化

上のフローチャートには書いていませんが、遺伝的アルゴリズムを実行して
実際にどのような経路なのか確認するためにmatplotlibを使って可視化します。

まずは都市の配置の可視化から。

実装例: show_cities

def show_cities(cities):
    for i in range(len(cities)):
        plt.scatter(cities[i][0], cities[i][1])

使用例: show_cities

cities = generate_rand_cities(10)
show_cities(cities)

使用例: show_cities 実行結果

Fig. 13 show_citiesで10都市の配置を可視化した例

次はルートの可視化。

実装例: show_route

def show_route(cities, genes):
    for i in range(len(genes)-1):
        if i == 0:
            plt.text(cities[int(genes[i])][0], cities[int(genes[i])][1], "start")
        else:
            plt.text(cities[int(genes[i])][0], cities[int(genes[i])][1], str(i))
        plt.plot([cities[int(genes[i])][0], cities[int(genes[i+1])][0]], 
                 [cities[int(genes[i])][1], cities[int(genes[i+1])][1]])
    plt.text(cities[int(genes[i+1])][0], cities[int(genes[i+1])][1], "goal")

使用例: show_route

cities = generate_rand_cities(10)
genes = generate_init_genes(10, 10)
show_route(cities, genes[0])

使用例: show_route 実行結果

Fig. 14 show_routeで巡回する順番を可視化

統合

これまでの関数を使って遺伝的アルゴリズムでTSPをとくプログラムです。

パラメータ

num_cities = 20
indivisuals = 21
generation = 10000
elite = 9
p_mutation = 0.005

Table 1 TSP_GAプログラムにおけるパラメータ設定

項目変数名例での値
都市数num_cities20
個体数indivisuals21
世代数generation10000
エリート選択数elite9
突然変異率p_mutation0.005

初期化

cities = generate_rand_cities(num_cities)
genes = generate_init_genes(indivisuals, num_cities)
show_cities(cities)

Fig. 13のような形で初期設定の都市の配置が確認できます。

実行部

top_indivisual=[]
max_fit = 0
for i in range(generation):
    fitness_vec = np.reciprocal(genes_path(genes, cities))
    child = np.zeros(np.shape(genes))
    for j in range(int((indivisuals-elite)/2)):
        parents_indices = roulette_choice(fitness_vec)
        child[2*j], child[2*j+1] = partial_crossover(genes[parents_indices[0]], 
                                                     genes[parents_indices[1]])
    
    for j in range(indivisuals-elite, indivisuals):
        child[j] = genes[np.argsort(fitness_vec)[j]]

    child = translocation_mutation(child, indivisuals-elite, p_mutation)
    top_indivisual.append(max(fitness_vec))
    genes = child
    if max(fitness_vec) > max_fit:
        max_fit = max(fitness_vec)
        show_cities(cities)
        show_route(cities, child[np.argmax(fitness_vec)])
        plt.text(0.05, 0.0, "generation: "+str(i))
        plt.show()

上で紹介していない要素について説明を追加します。
1行目のtop_indivisualはその世代のトップの適応度を格納していきます。
そうすると、以下のようなコードで世代に対しての適応度の変化を
確認することができます(Fig. 15)。

plt.plot(top_indivisual)
Fig. 15 世代ごとのトップ個体の適応度

11、12行目はエリート選択を行っています。
numpy.argsortで適応度が高い個体のインデックスを拾ってきています。

17~22行目まではわかりやすいように、現世代のトップ個体が
ハイスコアを出した場合にルートを可視化しています(Fig. 16)。

Fig. 16 実行部で出力されるルートを可視化した画像(抜粋)

5. まとめ

交叉の方法や突然変異の方法を別のに変えて、
またやってみたいと思いました まる

プログラム(GA_TSP.ipynb)公開リンク
Google Colab
Github


ロボロフスキーハムスター

珍しく多頭飼い可能なハムスター

小柄で人に懐きにくい性格

ハムスターカフェmogumoguにて撮影

「pythonで遺伝的アルゴリズム(GA)を実装して巡回セールスマン問題(TSP)をとく」への1件のフィードバック

  1. ピンバック: モンテカルロ法を使ってみる(円周率の計算) | 有閑是宝

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です