巡回セールスマン問題をいろんな近似解法で解く（その5: 焼きなまし法）

# はじめに

本記事は、[巡回セールスマン問題をいろんな近似解法で解く（その4: 反復局所探索法）](https://qiita.com/take314/items/8cdaf32e761665a54799)の続きとなります。

### 他の記事 （2022/05/30 追記）
他の記事へのリンクを以下に記載しましたので、ご興味ある方はぜひ読んでみてください。

- [その1: 総当たり法とグリーディ法](https://qiita.com/take314/items/dc2e6cf6d97889923c8b)
- [その2: 局所探索法と2-opt近傍](https://qiita.com/take314/items/33843f980c4cbab140ac)
- [その3: ランダム化最近近傍法とGRASP法](https://qiita.com/take314/items/cf12c29aef7bf8e43eeb)
- [その4: 反復局所探索法](https://qiita.com/take314/items/8cdaf32e761665a54799)
- [その6: さまざまな近傍](https://qiita.com/take314/items/69b93481403feb857d6e)

# 焼きなまし法

これまで紹介してきた手法は、基本的には近傍解の評価時に改善解への移動のみを許すものでした。今回紹介する**焼きなまし法**では、探索が進むにつれて下がっていく**温度**と呼ばれるパラメータを導入し、温度と解の改悪幅に応じた確率で改悪解への移動も許容することで、探索の序盤では多様な探索を、終盤では集中的な探索を行うことができます。

## 概念図

以下に概念図を示します。
<img width=90% src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/283041/1b17ffb4-d6f1-12c4-6d2a-4d0a161fa955.jpeg">

横軸は解の候補のバリエーションを表し、縦軸は総距離（TSPの解の評価値）を表します。実際の解空間は高次元すぎて図示することが不可能ですが、この図では簡単のため2次元で表しています。グラフを見ると山や谷があることがわかりますが、谷はそれぞれ局所最適解（そのうち最も低いものは大域最適解）に対応しています。また、ある谷から別の谷に移動するためにはその間にまたがる山を越えなければなりません。このような移動を行うためには一時的に総距離が増える方向、つまり改悪解への移動を許容する必要があります。

焼きなまし法では、温度の高い序盤の探索（暖色系の点線矢印）においては改悪解への移動が起きやすい、つまり大きな山を乗り越えやすいため、局所最適解にとらわれることなく解空間を幅広く探索できます。また、温度の低い終盤の探索（寒色系の点線矢印）では山を超える力は小さくなっていくものの、逆に言うと谷底の周辺を綿密に探索する力が大きくなっていくため、最終的にはある局所最適解かそれより僅かに悪い程度の解に到達することが可能となります。さらに、理論的には温度の下げ方を上手く選ぶと大域最適解に到達可能であることが保証されています。


## 擬似コード

次に、基礎的な焼きなまし法の疑似コードを以下に示します。

```simulated_annealing.py
# 焼きなまし法のpython擬似コード
def simulated_annealing(x_0 = 初期解, t_0 = 初期温度, c = 冷却率):
    x_current = x_best = x_0  # x_current, x_bestを初期解で初期化
    t = t_0  # 温度を初期化

    while 終了条件を満たさない:
        x = x_currentの近傍からランダムに選んだ解
        d = f(x) - f(x_current)  # f: 評価関数
        if d <= 0:
            x_current = x  # 改善解ならば常に移動する
            if f(x) < f(x_best):
                x_best = x
        elif rand() <= e ** (-d / t):  # e: ネイピア数, rand: 0以上1以下のランダムな実数を返す関数
            x_current = x  # 改悪解であっても、温度と改悪幅に応じた確率で移動する
        t = c * t  # 冷却する

    return x_best
```

$d$は解の**改悪幅**で、移動後の解の評価値から移動前の解の評価値を引いた値になります。いまはルートの総距離の最小化問題を考えており値が小さいほど良い解となるため、$d$が正の場合は移動後の解のほうが評価値が大きい、つまり解の改悪を表しますが、負の場合は逆に改善を表すことになります。

$t$は温度であり正の実数です。$c$は**冷却率**と呼ばれるパラメータで、0より大きく1未満の実数である必要があります。ループが回る毎に温度$t=ct$として更新されることから、探索が進むにつれて温度を下げる効果があることがわかります。ちなみにこのような温度の下げ方を**幾何冷却法**と呼びます。

改悪解への移動確率を表す指数関数$e^{-d/t}$については、指数部分の符号がマイナスであることから、温度$t$が大きいほど、あるいは改悪幅$d$が小さいほど指数関数全体としては大きくなり、その値は1に漸近します。このような性質により、焼きなまし法では探索の序盤においては温度が高く、改悪解への移動確率が大きくなるため、ある解領域から別の遠く離れた解領域への移動が起きやすく、探索の多様化力が強いといえます。逆に温度が低くなる探索の終盤では改悪解に移動しにくくなるため、局所最適解周辺の解を綿密に探索するようになり、集中化力が強いといえます。

ちなみに$t→∞$の極限においては改悪解への移動が常に許されるため、常にランダムな近傍解へと移動し続けることになります。一方で$t→0$の極限では改悪解への移動は常に却下されるようになるため、（近傍解の評価順序の差異を除き）通常の局所探索法と同様の動作をするようになります。

また細かな補足ですが、指数関数の底として必ずしもネイピア数を用いる必要はなく、必要であれば任意の底への変換を行い、そこに生じた差分を温度を定数倍して吸収することで等価な式を作れますので、それで計算を行ってもかまいません。軽く調べたところ、ネイピア数を使う流儀は物理学の文脈におけるメトロポリス法と呼ばれる計算手法が起源のようです。

そのほかの注意点として、今回は初期温度や冷却率をアルゴリズムの入力として定数で与えていますが、探索の初期段階で適切な値を推定する手法もあります。また理論的に大域最適解への到達を保証する**対数冷却法**や、ランダムな近傍解を使わない手法、終了条件として「十分に温度が下がった」ことを検知する手法もあったりなど、基本となる構造がシンプルであるがゆえに様々な変形や拡張が考案されています。詳しくは参考文献をご参照ください。

## 実験

それでは近傍に2-optを組み込んだ焼きなまし法をいくつかの問題例に適用してみます。使用した計算機はMacBook Pro（1.4 GHz クアッドコア Intel Core i5, 16 GB RAM）です。

...

### 48都市TSP
この問題例は最適値（すべての解の中で最も小さい総距離）が求められており、その総距離は10628です[^1]。
<img width=80% src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/283041/514ebf8f-1789-6948-535c-49e7e065ef3e.png">

総移動距離は10711であり、最適値に対しておよそ+0.78%の結果。計算時間は0.012秒であった。

...

###  1379都市TSP
この問題例も最適値が求められており、その総距離は56638です[^1]。
<img width=80% src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/283041/0428dee1-f139-5415-39f7-2c8b1bc4da05.png">

総移動距離は58572であり、最適値に対しておよそ+3.4%という結果。計算時間は約4.0秒であった。これまでの連載でのベスト値であったGRASPの60512を大きく更新できた。


## これまでの手法の結果と比較

### 表
|  |att48(10628) |kroA100<br>(21282) |tsp225<br>(3916) |att532<br>(27686) |nrw1379<br>(56638) |
|---|---|---|---|---|---|
|その1: 最近近傍法 |12861<br>0.000秒 |27807<br>0.000秒 |5030<br>0.000秒 |35516<br>0.000秒 |68964<br>0.000秒 |
|その2: 局所探索法<br>（最近近傍法+2-opt）  |11010<br>0.002秒 |22399<br>0.002秒 |4304<br>0.013秒 |29739<br>0.192秒 |61061<br>0.435秒 |
|その3: GRASP法<br>（ランダム化最近近傍法+2-opt）  |10895<br>0.013秒 |**21843**<br>0.014秒 |4158<br>0.103秒 |29320<br>1.142秒 |60512<br>5.009秒 |
|その4: 反復局所探索法<br>（2-opt+double bridge）  |10739<br>0.013秒 |22177<br>0.022秒 |4208<br>0.099秒 |29718<br>1.036秒 |61061<br>4.317秒 |
|今回: 焼きなまし法<br>（2-opt）  |**10711**<br>0.012秒 |22190<br>0.016秒 |**3988**<br>0.060秒 |**28591**<br>0.834秒 |**58572**<br>3.966秒 |

これまでの連載で紹介してきた手法の結果をまとめた表です。問題例名の下の丸括弧は最適値を表しています。また、各問題例においてこれまで紹介した手法の中でベストな総距離の値を太字で示しています。さらに、今回から総距離の下に計算時間も記載しました（すべて計算し直したため、過去の結果と多少の誤差が生じています…）。

### 結果と考察
今回は主に過去記事で紹介したGRASP法と反復局所探索法との比較を行いたかったため、それらのアルゴリズムにかかった秒数と同程度（か多少短い程度）となるようにそれぞれの問題例で終了条件と冷却率を調整した上で実験しました。かなりざっくりとした比較実験のため参考程度に留めておくべき結果だとは思いますが、5つの内4つの問題例でベストを更新できていることがわかります。

また、なんとなくではありますが都市数が多い問題例ほど更新幅が大きい傾向があるようにも見えます。これについては、過去記事においてGRASP法と反復局所探索法は終了条件を10イテレーション固定としていたため、いずれの問題例でも10個の局所最適解しか評価できておらず、大規模な問題例になるにつれて解空間の広さに対する探索の多様性が足りなくなっていたのかもしれません。それに対して今回の焼きなまし法では、探索中に局所最適解に到達する保証がないものの、近傍解の評価のたびに改悪解への移動が起きる可能性があるため、前者2つの手法と比べて探索の多様性が格段に大きくなったことが効いているのかな、と感じます。


# まとめ
焼きなまし法においては、温度と解の改悪幅に応じた確率で改悪解への移動を許容することで、探索の序盤では多様化を、探索の終盤では集中化を強めることがわかりました。またそれにより、比較的シンプルな実装にもかかわらず良質な解を得られることがわかりました。

次回はこれまでの連載で紹介できなかった2-opt以外の近傍や、それを用いた探索の手法をご紹介できればと思います（最終回の予定です）。


# 参考文献
柳浦睦憲 / 茨木俊秀 著、『[組合せ最適化 -メタ戦略を中心として-](http://www.asakura.co.jp/books/isbn/978-4-254-27512-4/)』

[^1]: TSPLIB: https://www.iwr.uni-heidelberg.de/groups/comopt/software/TSPLIB95/


# おわりに
ここまで読んでいただきありがとうございました。
株式会社オプティマインドでは、一緒に働く仲間を大募集中です。
カジュアル面談も大歓迎ですので、気軽にお声がけください。

【エンジニア領域の募集職種】
●[ソフトウェアエンジニア](https://hrmos.co/pages/optimind/jobs/0002012)
●[QAエンジニア](https://hrmos.co/pages/optimind/jobs/0002011)
●[Androidアプリエンジニア](https://hrmos.co/pages/optimind/jobs/0002008)
●[組合せ最適化アルゴリズムエンジニア](https://hrmos.co/pages/optimind/jobs/0002013)
●[経路探索アルゴリズムエンジニア](https://hrmos.co/pages/optimind/jobs/0002014)
●[バックエンドエンジニア](https://hrmos.co/pages/optimind/jobs/0002015)
●[インフラエンジニア](https://hrmos.co/pages/optimind/jobs/0002016)
●[UXUIデザイナー](https://hrmos.co/pages/optimind/jobs/0002010)

【ビジネス領域の募集職種】
●[セールスコンサルタント](https://hrmos.co/pages/optimind/jobs/0002001)
●[オペレーションコンサルタント](https://hrmos.co/pages/optimind/jobs/0002027)



『オプティマインドってどんな会社？』については、[こちら](https://www.optimind.tech/recruit/)から
Wantedlyでも[こちらで募集中](https://www.wantedly.com/companies/optimind/projects)