「道の駅スタンプラリー」をアルゴリズムで解く

## 問題提起

エンジニアのみなさん、道の駅に行ったことはありますか？
以後あるものと仮定します。

最近では、単にいつでもトイレが使えて車を止めて休憩できる施設(登録の必須要件です！)にとどまらず、直売所やレストラン、温泉、アクティビティ施設などを併設した魅力的な道の駅が増えています。

北海道には、現在127の道の駅が登録されており、毎年度、その道の駅を巡るスタンプラリーが開催されています。しかし、北海道全土に点在する127もの道の駅を1年で巡るのはなかなかに大変です。何かいい方法はないでしょうか？

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/a6bedc84-f2ef-3a3c-3548-da65957c762b.png)


すでに多くの人が気づいているかもしれませんが、これは有名な巡回セールスマン問題と同じ格好をしています。巡回セールスマン問題は、その厳密な最適解を求めることは現時点では不可能に近いとされていますが、比較的よい解を求める方法は色々と考案されています。

今回は、巡回セールスマン問題の近似解法を使って、道の駅スタンプラリーを最適化してみたいと思います。
記事の前半では、最適化の前段階として、データの取得や表示部分の実装を解説します。
後半では、巡回セールスマン問題の解法と道の駅スタンプラリーにおける実例を解説します。

## データの取得

競技プログラミングであれば、回るべき点同士の距離関係は入力として与えられますが、今回は自力で取得する必要があります。

実は、今回の道の駅のような、日本の地理情報や公共施設に関係する統計データは、**オープンデータ** として、国や自治体のウェブブラウザ上で機械可読な形式で積極的に公開されています。

ということで、道の駅に関するオープンデータが存在するかどうか検索してみたところ、やはり存在したのでこれを活用してみたいと思います。

https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-P35.html

便利なことに、上記データから北海道に限定したデータも配信されていたため、そちらを使うことにします。

https://www.geospatial.jp/ckan/dataset/ksj-p35-01

このデータは geojson 形式で構造化されており、つまるところ、特定のスキーマにそったただのjsonなので普通にparseできます。

と思ったのですが、データを見てみると全部で122件しかデータがないことがわかります。どうやら最近更新されていないようで、ここ数年の新規開設と廃止の動きが反映されていないようなので、差分データを手入力で作成することにしました。

余談ですが、今リンクを取得しに github を訪問したところ、地図上でgeojsonがプレビューできていました。そんな機能あったんですね。

https://github.com/2427dkusiro/MichinoekiTSP/blob/master/MichinoekiTSPDataLib/Resources/Michinoeki_diff.json

https://github.com/2427dkusiro/MichinoekiTSP/blob/master/MichinoekiTSPDataLib/Resources/Michinoeki_Abolish_List.txt

これで、127の道の駅に関する基本情報は収集できました。

## 経路の取得

### 方針

さて、道の駅スタンプラリーを最適化するためには、道の駅間の移動コストを計算する必要があります。今回は、道路や交通状況を加味した所要時間を google maps platform を活用して計算することにします。
なお、こんなことに通行料を払うのはアホくさいので、有料道路は使用しない経路を探索します。

今回は、最新の Routes API を使用して計算することにします。

https://developers.google.com/maps/documentation/routes?hl=ja

こちらのAPIは、リクエスト内容に応じて1回あたり 0.005ドル ~ 0.015ドル の料金設定となっています。比較的高額な請求となる可能性があるため、注意して使用してください。なお、google maps apiでは、月200ドルの無料枠と、初回登録時にgoogle cloudのサービス全体で使用できる300ドルの無料枠が設定されています。

今回のプロジェクトでは、127の道の駅同士の組み合わせのため、127x126 = 16002 件のリクエストになることが想定されます。これは、失敗しなければ無料枠で収まる計算です。

ほとんどダウンロードされていないのが不安になりますが、今回コーディングに使用するC#用のRoutes APIクライアントが公式から提供されていますので、これを使用します。

https://www.nuget.org/packages/Google.Maps.Routing.V2/1.0.0-beta09

### 実装

どこにも資料がないので軽く使い方を説明しておきます。興味がない人は次のセクションまで飛ばしてください。
コード全文は以下にあります。

https://github.com/2427dkusiro/MichinoekiTSP/blob/master/MichinoekiTSPDataLib/GoogleRouteAPI.cs

まず、認証情報を環境変数に設定します。google maps platformの管理コンソール上で新しいプロジェクトを作成し、`google_credential.json` を適当な位置にダウンロードしておきます。今回は、C#プロジェクトのルートに配置しました。

環境変数 `GOOGLE_APPLICATION_CREDENTIALS` に、先程のファイルのパスを設定します。入りっぱなしになっているのは嫌な気がするので、コンストラクタで都度設定することにします。

```cs
public GoogleRouteAPIClient()
{
    var execPath = Path.GetDirectoryName(Assembly.GetExecutingAssembly().Location);
    var credential = Path.Combine(execPath!, "google_credential.json");
    Environment.SetEnvironmentVariable("GOOGLE_APPLICATION_CREDENTIALS", credential);
}
```

次に、経路を取得する部分です。取得を希望する項目を `X-Goog-FieldMask` ヘッダーを使用して指定します。今回はひとまず全取得にしました。

次に、リクエスト本文を作成します。
出発地 `Origin` と到着地 `Destination` の緯度経度を指定し、使用する交通機関 `TravelMode` を指定します。今回は自動車です。
`RoutingPreference` はリアルタイム交通状況を加味するかどうかです。API使用料金に影響するようです。
`PolylineQuality` は経路を表す線 `Polyline` の解像度です。解像度が高い方が曲線がなめらかになり実際の道路形状とルート形状が一致しますが、データ量が増えます。ひとまず高解像度にしましたが、この程度の用途であれば明らかに過剰でした。
`RouteModefiers` はおなじみのgoogleマップでチェックボックスで指定するやつです。今回は、有料道路とフェリーを禁止して、高速道路は許可します。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/987790b0-7951-08c7-e6b0-3b84d2324599.png)

コードは以下のようになります。

```cs
using Google.Api.Gax.Grpc;
using Google.Maps.Routing.V2;
using Google.Type;

RoutesClient _client = await RoutesClient.CreateAsync();
var callSettings = CallSettings.FromHeader("X-Goog-FieldMask", "*");
var request = new ComputeRoutesRequest
{
    Origin = new Waypoint
    {
        Location = new Location { LatLng = new LatLng { Latitude = from.Latitude, Longitude = from.Longitude } }
    },
    Destination = new Waypoint
    {
        Location = new Location { LatLng = new LatLng { Latitude = to.Latitude, Longitude = to.Longitude } }
    },
    TravelMode = RouteTravelMode.Drive,
    RoutingPreference = RoutingPreference.TrafficUnaware,
    PolylineQuality = PolylineQuality.HighQuality,
    RouteModifiers = new RouteModifiers()
    {
        AvoidFerries = true,
        AvoidTolls = true
    }
};
```

次に、実際にリクエストを投げます。しかし、時々APIコールが失敗するため、リトライ処理を書いています。
リトライするときは、ある程度待ってからクライアントを再作成して投げます。

```cs
ComputeRoutesResponse? response = null;
int retry = 0;

async Task TryGet()
{
    try
    {
        response = await _client.ComputeRoutesAsync(request, callSettings) ?? throw new InvalidOperationException("API response was null.");
        if (!response.Routes.Any())
        {
            throw new InvalidOperationException("API returned an invalid response.");
        }
    }
    catch (Exception ex)
    {
        writeLog?.Invoke($"API Error Response:{ex}");
        if (retry++ < maxRetry)
        {
            await Task.Delay(1000);
            _client = await RoutesClient.CreateAsync();
            await TryGet();
        }
        else
        {
            throw;
        }
    }
}
await TryGet();
```

次に結果を解析します。いつもの googleマップ 同様、複数の候補ルートが返却できるようなレスポンス構造になっていますが、今回は最適ルートのみを使用します。再度になりますが、`Polyline` というのは経路の形を表す線です(通る座標点の集合です)。

```cs
Google.Maps.Routing.V2.Route route = response.Routes.First();
var title = route.Description;
var distance = route.DistanceMeters;
var duration = TimeSpan.FromSeconds(route.Duration.Seconds + (route.Duration.Nanos / 1_000_000_000));
var averageSpeed = distance / duration.TotalHours * 1000;
var polyline = route.Polyline.EncodedPolyline;
```

上記コードで解析したレスポンスの一例です。(レスポンス原文ではありません)

```json
{
"From":"☆ロマン街道しょさんべつ",
"To":"230ルスツ",
"Title":"国道232号 と 国道275号",
"DistanceMeters":281617,
"Duration":"05:31:38",
"Polyline":"{w~nGoxi`ZCOGy@COC[Co@?c@B]@a@FcAD(以下略)"
}
```

このパッケージでは、HTTP3による呼び出しができるようで、またgoogle maps apiは結構高速なこともあり、一回1円以上するAPIですがかなりのスピードで呼び出しができます。1時間程度ですべての道の駅ペア間の計算が完了しました。

### 後日談

今回、交通状況を加味しない (`TrafficUnaware`) を指定したので、Routes API のドキュメントを読む限りでは、0.005ドル / req になるはずです。

https://developers.google.com/maps/documentation/routes/usage-and-billing?hl=ja

しかし、後日管理コンソールで確認すると、0.015ドル / req で請求されていました。怖いです。
予想外にも、google maps api の無料枠をはみ出しています。google cloud全体の初回無料枠で踏みとどまったのでよかったですが、かなりホラー案件です。~~やっぱりクラウドこわい~~

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/f2906133-b20f-05b4-927c-2fb4e1d19c4f.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/35cad08d-2c24-1ea3-edf9-4d652edc0634.png)

## ビジュアライザをつくる

### UI実装

さて、ここから経路探索といきたいところなのですが、先に探索した経路を可視化するビジュアライザをつくりましょう。やはり、地図上で表示されているほうがわかりやすいですからね。

先程の Routes API のコールや、このあとの経路探索部の実装はネイティブでやるのですが、地図表示部分だけはwebでやります。アプリケーションとしては、windowsのGUIアプリケーションとして、MicrosoftのWebview2で表示を埋め込みます。ネイティブ側からは、 `ExecuteScriptAsync` 関数で Top-level のJS関数を呼び出せますので、相互の連携も可能です。

web側は、OpenStreetMap の地図データを使用し、Leaflet というライブラリを使用して表示します。
地図上にマーカーや経路の線を表示する機能がライブラリから提供されているため、これをネイティブ側から呼び出すように明示的にラップして、top-levelにおいておきます。正直、コードを見てもらったほうが早いと思うので貼っておきます。

```html
<!DOCTYPE html>
<html>

<head>
    <title>地図表示</title>
    <meta charset="utf-8">
    <meta name="viewport" content="width=device-width, initial-scale=1.0">
</head>

<body>
    <link rel="stylesheet" href="https://unpkg.com/leaflet@1.6.0/dist/leaflet.css"
          integrity="sha512-xwE/Az9zrjBIphAcBb3F6JVqxf46+CDLwfLMHloNu6KEQCAWi6HcDUbeOfBIptF7tcCzusKFjFw2yuvEpDL9wQ=="
          crossorigin="" />
    <!-- Make sure you put this AFTER Leaflet's CSS -->
    <script src="https://unpkg.com/leaflet@1.6.0/dist/leaflet.js"
            integrity="sha512-gZwIG9x3wUXg2hdXF6+rVkLF/0Vi9U8D2Ntg4Ga5I5BZpVkVxlJWbSQtXPSiUTtC0TjtGOmxa1AJPuV0CPthew=="
            crossorigin=""></script>

    <div id="main"></div>

    <style>
        * {
            margin: 0px;
            padding: 0px;
        }

        #main {
            width: 100vw;
            height: 100vh;
        }
    </style>

    <script lang="javascript">
        function delay(time) {
            return new Promise(resolve => setTimeout(resolve, time));
        }

        globalThis["ensureInitialized"] = async () => {
            while (globalThis["mapObj"] === undefined) {
                await delay(100);
            }
        }

        const map = L.map('main');//.setView([35.66572, 139.73100], 15);
        const tileLayer = L.tileLayer('https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
            attribution: '© <a href="http://osm.org/copyright">OpenStreetMap</a> contributors, <a href="http://creativecommons.org/licenses/by-sa/2.0/">CC-BY-SA</a>',
        });
        tileLayer.addTo(map);
        globalThis["mapObj"] = map;

        function setInterop() {
            globalThis["setView"] = (lat, lng) => {
                globalThis["mapObj"].setView([lat, lng]);
            }

            globalThis["setZoom"] = (level) => {
                globalThis["mapObj"].setZoom(level);
            }

            globalThis["addMarker"] = (lat, lng, name) => {
                const marker = L.marker([lat, lng]);
                marker.bindPopup(name);
                marker.addTo(globalThis["mapObj"]);
            };

            globalThis["addPolyline"] = (position, colorCode) => {
                if (globalThis["polylines"] === undefined) {
                    globalThis["polylines"] = [];
                }
                const line = L.polyline(position, {
                    "color": colorCode,
                    "weight": 5,
                    "opacity": 0.6
                }).addTo(globalThis["mapObj"]);
                globalThis["polylines"].push(line);
                return globalThis["polylines"].length - 1;
            }

            globalThis["removePolyline"] = (id) => {
                const line = globalThis["polylines"][id];
                globalThis["mapObj"].removeLayer(line);
            }
        }
        setInterop();
    </script>
</body>
</html>
```

### Polylineの相互運用

ただし、1点注意事項があります。google のいうpolylineと、leafletのいうpolylineのエンコード形式が一致しません。というのも、google から取得したpolylineは独自にエンコードされているからです。

すでに、この形式のencoder/decoderの実装が存在したため、今回はそれを使用しています。先人に感謝です。

https://gist.github.com/shinyzhu/4617989

なお、オリジナル版の `CoordinateEntity` 型が何を指しているのか判別できなかったため(少なくとも標準で `using` が通る範囲にはない) ため、少しいじっています。

適当な `record` 型を用意して、それに差し替えました。

```cs
/// <summary>
/// 地図上の点を表現します。
/// </summary>
/// <param name="Name">この地点の名前。</param>
/// <param name="Latitude">緯度。</param>
/// <param name="Longitude">経度。</param>
public record GeometryPoint(string? Name, double Latitude, double Longitude);
```

変更後のコードです。

https://github.com/2427dkusiro/MichinoekiTSP/blob/master/MichinoekiTSPDataLib/PolylineEncoder.cs

あとは、JS文法の配列な文字列に変換してJS側に渡せば経路が表示されます。

```cs
public async Task AddPolyline(IEnumerable<GeometryPoint> points, string color)
{
    var arrayString = $"""[{string.Join(',', points.Select(x => $"[{x.Latitude},{x.Longitude}]"))}]""";
    var result = await webView.CoreWebView2.ExecuteScriptAsync($"""addPolyline({arrayString},"{color}")""");
    var id = JsonSerializer.Deserialize<int>(result);
    polylines.AddLast(new Polyline(id, points));
}
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/596167c8-2ab8-3997-d5c9-6b675d58ebaf.png)

## 経路探索アルゴリズム

### 概要

さて、道の駅を巡る効率のいい経路を計算するにはどうすればいいでしょうか。道の駅を巡る経路とは、言い換えればどの順番で道の駅を訪問するかということになります。すべての訪問する順番の選択肢に対して、所要時間を合計し、一番少ないものを選択すればそれは間違いなく最適解なのですが、訪問する順番の選択肢とは127ある道の駅を並べる順列であり、`127!` 通りあります。これは、10の213乗程度であり、とても計算できないほどに膨大です。

現時点で、この問題を解決する方法は見つかっていません。若干マシな方法はあるようなのですが、若干マシ程度であって、やはり計算できないほど膨大な組み合わせに爆発します。ちゃんと説明するならば、多項式時間でこの問題を解くアルゴリズムは見つかっておらず、また見つかる望みもほとんどないというわけです。

そこで、「最善」かどうかはわからないが、限りなく最善に近い経路を見つける方法を考えます。現実的にはかなりいい妥協点なわけです。こうした比較的よい解を探す方法は様々なものが考案されており、巡回セールスマン問題(TSP)の解法をキーワードに調べれば、様々なものを知ることができるでしょう。

アルゴリズムに関しては、以下の記事などを参考にしながら、

https://future-architect.github.io/articles/20211201a/#%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0

- Nearest Neighbor
- 2-opt LS (実装上の工夫なし)
- 2-opt ILS (実装上の工夫なし)
- 2-opt SA

を実装します。雰囲気を掴んでもらえるようにコード断片を貼っていきますが、道の駅クラスや経路クラスなど断り無く使用しますので、全体は前述のリポジトリから確認してください。なお、経路探索では、道の駅登録番号が最も若い１番で、札幌からもほど近い「三笠」の道の駅を起点に、単純な運転時間の総和が最も少なくなるように最適化します(道の駅に立ち寄る時間や営業時間は考慮しない)。

### Nearest Neighbor

#### 概要と実装

Nearest Neighbor法は、貪欲に未到達の道の駅の中から一番近い道の駅に向かうことを繰り返す方法です。そのため、実装は非常に簡単です。しかし、最初はきれいな経路が描かれますが、後のことを考えずに道の駅を巡っていくので、最後の方の経路が汚くなることが予想されます。

```cs
public TSPAnswer Solve()
{
    Dictionary<GeometryPoint, bool> stationDict = context.Michinoekis.ToDictionary(x => x, _ => false);
    int count = stationDict.Count;
    if (stationDict.ContainsKey(context.StartPoint))
    {
        stationDict[context.StartPoint] = true;
        count--;
    }
    var ans = new Route[count];

    static Route NonVisitedMin(Route[] routes, Dictionary<GeometryPoint, bool> stationDict)
    {
        Route? min = null;
        foreach (var r in routes)
        {
            if (stationDict[r.To])
            {
                continue;
            }
            if (min is null)
            {
                min = r;
                continue;
            }
            if (min.Duration > r.Duration)
            {
                min = r;
                continue;
            }
        }
        return min!;
    }

    var current = context.StartPoint;
    for (int i = 0; i < ans.Length; i++)
    {
        var selected = NonVisitedMin(context.FindEdgeFrom(current), stationDict);
        current = selected!.To;
        stationDict[current] = true;
        ans[i] = selected;
    }

    return context.SubmitAnswer(ans);
}
```

#### 結果

結果を示します。

| アルゴリズム | 所要時間 | 距離 |
|:-:|:-:|:-:|
| Nearest Neighbor |  3日14時間53分08秒 | 4671km+304m |

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/8ab61678-c9d9-27b2-a715-22824aa58ba8.png)

解を分析してみましょう。
まず道東、網走付近の経路を拡大します。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/fff79bc5-c36b-2665-724a-5f02fce06661.png)

最短を繰り返した結果、行ったり戻ったりしてしまい、経路がごちゃごちゃしていることがわかります。
次に、日高付近を拡大します。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/643e1806-d0df-9c1b-1d5b-3de90d7a2c53.png)

線が重なるのでわかりにくいですが、完全に同じ道を戻っています。
といった具合に、最初はよくてもそれが後から困るケースが出てくるため、Nearest Neighborではあまり効率的な経路が描けてはいません。

### 2-opt LS

#### 2-optとは？

Nearest Neighborではあまりよい解が得られないことがわかりました。
ここで、この完璧ではない解を修正・改良して、よりよい解を得ることを考えてみましょう。
改良するには、道の駅を訪問する順番を入れ替えればよいわけです。

厳密解を得るには、127個の道の駅を互いに入れ替える方法を考える必要がありますが、それが不可能なのは説明した通りですから、数個の道の駅を選択して、それを所要時間が短くなるように入れ替える操作を考えてみましょう。この方法では、厳密解を得ることはできませんが、2個の道の駅を選択して所要時間が短くなるように入れ替えることを繰り返す場合、バブルソートが `O(N^2)` であるように、この操作は `N=127` において現実的な計算量に収束します。この手法を、2個を選択して最適化することから **2-opt** と呼びます。

#### 局所解

順列に対して、2点スワップを繰り返すという発想は同じですが、バブルソートでは厳密解を得ることができるのに対して、巡回セールスマン問題では厳密解を得ることができません。
実は、すべての道の駅が同じただ１本の道路上にあれば、その最適解はソートと一致します。1次元であれば、厳密解を得ることができるのです。しかし、現実には、`(緯度 , 経度 , 車の流れ)` からなる空間上に道の駅は存在しており、次のような問題が発生します。

2つの道の駅の順番を入れ替えることを繰り返すと、やがてこれ以上改良できない地点に到達します。しかし、この地点が最適とは限りません。なぜなら、一旦所要時間が伸びるような入れ替えをすることでしか見えない経路が存在する可能性があるからです。したがって、この方法では得た解は狭い範囲だけを見たときの最適解、局所解を得るというわけです。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/d7c5888e-1263-0e80-8968-33f9e2e04bcc.png)

※画像は借りました(https://zellij.hatenablog.com/entry/20120825/p1)

#### 実装

例えば、以下のような経路があるとします。

```
三笠
芦別
滝川
(略)
厚田
新篠津
当別
```

この `芦別` と `新篠津` を入れ替えます。このとき、`芦別 - 滝川 - (略) - 厚田 - 新篠津` はつながっているので、入れ替えるというよりは、この部分をすべてひっくり返すことになります。ひっくり返すとは、出発点と到着点をすべておのおの逆にするということです。

```
三笠
新篠津
厚田
(略)の反転
滝川
芦別
当別
```

元々の移動時間は、`三笠 -> 芦別` + `新篠津 -> 当別` + (略の部分) でした。
新しい移動時間は、`三笠 -> 新篠津` + `芦別 -> 当別` + (略の部分の反転) になります。
これらの所要時間をそれぞれ計算し、所要時間短縮になるのであればこの変更を受け入れます。
どうしてもごちゃっとしたコードになりがちですが、できるだけきれいにこの方針を実装します。

```cs
public static bool TrySwap(TSPSolverContext context, int i, int j, Span<Route> routes, Route[] swapTable)
{
    var (oldCost, newCost) = CalcSwap(context, i, j, routes, swapTable);

    if (oldCost > newCost)
    {
        CommitSwap(i, j, routes, swapTable);
        return true;
    }
    return false;
}

public static (TimeSpan oldCost, TimeSpan newCost) CalcSwap(TSPSolverContext context, int i, int j, ReadOnlySpan<Route> routes, Route[] swapTable)
{
    Debug.Assert(i != j);
    Debug.Assert(routes.Length == swapTable.Length);
    if (i > j)
    {
        (i, j) = (j, i);
    }
    Debug.Assert(j < routes.Length);

    var a = routes[i];
    var b = routes[j];
    var newRouteA = context.FindRoute(a.From, b.From);
    var newRouteB = context.FindRoute(a.To, b.To);

    swapTable[i] = newRouteA;
    var oldCost = a.Duration;
    var newCost = newRouteA.Duration;

    var routeSlice = routes[(i + 1)..j];
    var swapTableSlice = swapTable.AsSpan()[(i + 1)..j];
    for (int k = 0; k < routeSlice.Length; k++)
    {
        var target = routeSlice[^(k + 1)];
        swapTableSlice[k] = context.FindRoute(target.To, target.From);
        oldCost += routeSlice[k].Duration;
        newCost += swapTableSlice[k].Duration;
    }
    swapTable[j] = newRouteB;
    oldCost += b.Duration;
    newCost += newRouteB.Duration;

    return (oldCost, newCost);
}

public static void CommitSwap(int i, int j, Span<Route> routes, Route[] swapTable)
{
    if (i > j)
    {
        CommitSwap(j, i, routes, swapTable);
        return;
    }

    swapTable.AsSpan()[i..(j + 1)].CopyTo(routes[i..]);
}
```

あとは、バブルソートのように探索していきます。

```cs
for (step = 1; step <= p.MaxSteps; step++)
{
    bool hasSwap = false;
    for (int i = 0; i < routes.Length - 1; i++)
    {
        for (int j = i + 1; j < routes.Length - 1; j++)
        {
            hasSwap |= TSPUtil.TrySwap(context, i, j, routes, swapTable);
        }
    }
    if (!hasSwap)
    {
        context.WriteMessage($"2-opt LSは{step - 1}stepsで極大値に達しました");
        break;
    }
}
```

#### 結果

結果を示します。

| アルゴリズム | 所要時間 | 距離 |
|:-:|:-:|:-:|
| Nearest Neighbor |  3日14時間53分08秒 | 4671km+304m |
| 2-opt LS | 3日05時間33分52秒 | 4159km+756m |

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/e2de893b-8e59-c9e5-9b1c-aeb49230e9cf.png)

9時間以上改善することができました。先程拡大した網走付近を含む、道東部分を拡大してみましょう。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/aec3ea95-458f-a7e2-d2a8-27e931e10cae.png)

線同士の交差や重複が大きく減り、ルートがきれいになったことがわかります。わかりやすいように、ほぼ同じ位置のNearest Neighborでの解を掲出します。 

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/7aafe356-a677-cc95-5454-c6d1e668fe08.png)

しかし、2-opt LS でもまだ最適ではないことがわかるでしょうか。例えば、道東エリアではオホーツク海に沿って道の駅を回収しているときに内陸側の道の駅を回収していないため、太平洋側から大きく迂回して回収に来ることになっていますし、道北エリアでも国道40号線沿いの道の駅を回収していないためやはりあとから戻ってくることになっています。

といった具合に、局所解ではNearest Neighborで見られたような線の交わりや重複という露骨な無駄は解消することができるのですが、直前の考察のような、道の駅の存在範囲をグループ化・ライン化して、その後のルートをきれいにすることを考慮して必要ならば若干遠回りしながら回収していく、といった大局的戦略までは見られず、やはり「局所解」であるということがわかります。

### 2-opt ILS

#### 概要

これ以上解を改善するには、局所解から抜け出す必要があります。そこで、局所解に陥ったときに局所解を脱出する操作を加え、また新たな局所解を得、という操作を繰り返すことを考えます。

局所解を脱出する操作として、ランダムに4つの道の駅を選択し、その順番を入れ替える Double-Bridge を採用します。これは、ある程度整った2-opt LSの解を完全に破壊することなく(127のうち4点を入れ替えるだけなので)、従来の2-optではたどり着けなかった解の世界を見えるようにすることができます。

ここから、ランダム性が入ってきます。そのため、結果は試行ごとに若干変化します。
Double-Bridgeの実装例を示します。

```cs
public static Route[] Kick(TSPSolverContext context, Route[] routes, Random? random)
{
    Span<int> target = stackalloc int[4];
    CreateRandoms(target, routes.Length, random);
    target.Sort();

    var arr = new Route[4];
    arr[0] = context.FindRoute(routes[target[0]].From, routes[target[2]].To);
    arr[1] = context.FindRoute(routes[target[1]].From, routes[target[3]].To);
    arr[2] = context.FindRoute(routes[target[2]].From, routes[target[0]].To);
    arr[3] = context.FindRoute(routes[target[3]].From, routes[target[1]].To);

    int current = 0;
    var newArray = new Route[routes.Length];
    for (int i = 0; i < newArray.Length; i++)
    {
        if (target.Contains(current))
        {
            var obj = arr[target.IndexOf(current)];
            newArray[i] = obj;
            current = routes.Select((x, i) => (x, i)).First(x => x.x.To == obj.To).i + 1;
        }
        else
        {
            Debug.Assert((uint)current < routes.Length);
            newArray[i] = routes[current++];
        }
    }

    Debug.Assert(newArray.Select(x => x.From).Distinct().Count() == routes.Length);
    return newArray;
}
```

#### 結果

結果を示します。これ以降の解にはランダム性があるため、100回程度実行して最良の解を選択することである程度ランダム性を抑えています。この反復は独立しているため、マルチスレッド化して実行しています。

| アルゴリズム | 所要時間 | 距離 |
|:-:|:-:|:-:|
| Nearest Neighbor |  3日14時間53分08秒 | 4671km+304m |
| 2-opt LS | 3日05時間33分52秒 | 4159km+756m |
| 2-opt ILS | 3日03時間08分34秒 | 4010km+514m |

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/36580053-3bab-6ea6-6baf-61e8f4ee2dcd.png)

国道40号に沿った道の駅をうまく回収することで、道北部分のルートがスマートになっています。音威子府から中川や美深の道の駅への経路はほぼ折り返し経路ですが、これにより大局的に見ると一回来るだけで済むようになっており合理的です。折り返し経路をほぼ完全に消し去った2-opt LSに比べ、大局的な最適化のためにところどころ折り返し経路を許容するように変化していることがわかると思います。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/51b5bcaa-3a4c-97fc-7b78-091ebd09d13c.png)

### 焼きなまし法(2-opt SA)

#### 概要

焼きなまし(元ネタ)について知るより先に焼きなまし法について知っていたのですが、焼きなまし法は焼きなましに着想を得た最適化アルゴリズムであり、疑似アニーリング(SA)法とも呼ばれます。

現実世界の原子は、そのエネルギーのポテンシャルが低い方に一般に遷移します。エネルギーが低い方が安定だからです。そして、安定な物質は一般に丈夫で有用です。そこで、究極に安定な、規則正しい配列の結晶構造をした金属を得たい。しかし、ただの金属では「それなりに」安定な、それなりに規則正しくて多少乱雑な構造をしており、それなりに安定なためそこから結晶構造が変化することは自然にはありません。そう、局所解に嵌っているのです。

そこで、あえて金属を加熱し、エネルギーが高い状態にします。その金属に衝撃を与えると、粒子は高エネルギー状態にあるため自由に動き、構造を変化させることができます。少しずつ温度が減少していくため、粒子は移動しにくくなり、最終的にある程度安定な状態に落ち着きます。このとき、元の状態より安定なところに遷移する可能性があるのです。つまり、局所解を脱して大局的な最適解に遷移する見込みがあります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/d7c5888e-1263-0e80-8968-33f9e2e04bcc.png)

※再度ですが画像は借りました(https://zellij.hatenablog.com/entry/20120825/p1)

というのを、アルゴリズムで活用します。温度というパラメータが高い場合は、所要時間が改善しない2点道の駅スワップでも受け入れられる可能性が高まります。やがて温度は減少していき、温度の減少に合わせて結果が悪化するスワップが拒否される結果が高まります。そして、温度が常温になると、改善するスワップのみが受け入れられ、その局所解で停止します。

なお、この手法はかなり適用対象が広く、巡回セールスマン問題にとどまらず、様々な最適化問題に応用することができます。

#### 実装

実は、2-opt LSの時点でスワップの評価と適用を分けていたので、簡単に実装できます。
温度変化は線形としました。

```cs
public TSPAnswer Solve()
{
    var answer = initialSolver.Solve();
    Route[] routes = answer.Routes.ToArray();
    Route[] swapTable = new Route[routes.Length];

    var random = context.Random;
    double st = parameter.StartTemp;
    double et = parameter.EndTemp;
    // const double temp_factor = 0.999;
    Span<int> buf = stackalloc int[2];

    var maxIteration = parameter.MaxIteration;
    for (int i = 0; i < maxIteration; i++)
    {
        var rand = TSPUtil.CreateRandoms(buf, routes.Length, random);
        var (oldCost, newCost) = TSPUtil.CalcSwap(context, buf[0], buf[1], routes, swapTable);

        double temp = st + (et - st) * i / maxIteration; // Math.Pow(temp_factor, i);
        double prod = Math.Exp((oldCost - newCost).TotalSeconds / temp);
        var r = random.NextDouble();
        if (prod > r)
        {
            TSPUtil.CommitSwap(buf[0], buf[1], routes, swapTable);
        }
    }

    return context.SubmitAnswer(routes);
}
```

#### 結果

焼きなまし法では、温度などパラメータ調整の余地があります。
今回は、色々試しながら、開始温度8000度、終了温度10度、イテレーション2百万回を設定しました。
このアルゴリズムもランダム性があるため、試行を重ねています。

| アルゴリズム | 所要時間 | 距離 |
|:-:|:-:|:-:|
| Nearest Neighbor |  3日14時間53分08秒 | 4671km+304m |
| 2-opt LS | 3日05時間33分52秒 | 4159km+756m |
| 2-opt ILS | 3日03時間08分34秒 | 4010km+514m |
| 2-opt SA | 3日00時間49分02秒 | 3873km+368m |

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/9ccb1e24-093d-a8dd-f7bf-2094bf093344.png)

かなり所要時間が少なくなりました。拡大して個別に評価してみましょう。

道北部分

オロロンラインを北上する途中で、経路をそれて中川の道の駅を回収することで、その後の浜頓別～美深へのルートをきれいにしている。2-opt LSでは苦手だった箇所だが、日本海(オロロンライン)、国道40号、オホーツク海側と、3系統(奇数！)ある特徴をうまく掴んできれいに捌いていると言える。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/fb1e6252-0ab3-006c-bf49-7145cbbf4c53.png)

道東部分

起点の三笠付近は道の駅が密集しており、かたまりがつかみにくいが、神居古潭より東側の平野部を先にジグザクになっても回収しておくといういい判断ができている。また、道北から続いてオホーツク海に沿ったあと、サロマ湖で引き返して内陸側を回収するという判断、また国道452号というかたまりの捉え方は合理的。2-opt LSでみれらた網走付近の迂回や日高の長距離折返しも解消しており、とても整っている。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/fd254838-bb79-9a15-1b50-9166da35b96e.png)

道南部分

高速をうまく活用できる七飯の経路判断のほか、道の駅が密集する胆振地方でもうまく2つのラインで捌いており、合理的なルートであるといえる。終点を望羊中山とする判断も、札幌方面に道の駅がない孤立地帯であることを考えると適切だろう。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/143cc6b1-07ef-4e10-67e6-6ee9f266e89f.png)

と、3日1時間を切る、非常に優秀な解が得られました。経路を目視で評価しても冗長であると思われる箇所はなく、厳密解にかなり近いことが推測されます。

## 最後に

丸3日ちょっと運転し続ければ道の駅をコンプリートできることがわかりました。
しかし、そんなことは現実的には無理なので少しずつ分割しながら効率よくスタンプを集めることにしました。

焼きなまし法により、無駄のない経路が発見されたことで、どんどんスタンプが集まる...
そんな中、伊達の道の駅に到着しました。伊達の道の駅には、農産物直売所と歴史博物館が併設されています。
そして、歴史博物館に併設の工房では、焼きなまし(物理)を間近で見ることができます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510423/8cb0a2c5-7989-4afc-fdfa-41e3e5899eee.png)

コンピュータアルゴリズムと刀鍛冶、だいぶかけ離れた存在のようですが、こうしてつながることもあるのです。全く無用の知識などない、どんな過去でも記憶するに値する、と実感した気がします。