{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "remove-cell" ] }, "outputs": [], "source": [ "import sys\n", "\n", "sys.path.insert(1, \"../../\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 巡回セールスマン問題\n", "\n", "Amplify の使用例として、巡回セールスマン問題を Amplify を用いて解く方法を解説します。\n", "巡回セールスマン問題とは、いくつかの都市の集合が与えられたとき、ある都市から出発してすべての都市を 1 回ずつ訪れたあと最初の都市に戻ってくる巡回路のうち、最も長さが短いものを求める組合せ最適化問題です。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 巡回セールスマン問題の定式化\n", "\n", "以下では、都市の数を `NUM_CITIES` と表記します。\n", "\n", "### 変数\n", "\n", "まず、`NUM_CITIES` 個の都市のうちどの都市を訪れるかを表すために、 `NUM_CITIES` 個のバイナリ (0-1) 変数を用いることにします。都市 `i` を訪れる場合、`i` 番目のバイナリ変数を 1 、それ以外の変数を 0 とすることで表現します。\n", "\n", "たとえば都市数が 5 であり 5 つのバイナリ変数の値が `[0, 0, 0, 1, 0]` の場合、インデックス 3 の変数のみが 1 となっているので都市 3 を訪れることの表現となっています。このように、n 個の選択肢の中からどれを選ぶかを n 個のバイナリ変数のうち 1 つを 1 にすることで表す方法を one-hot エンコーディングといい、機械学習などでもよく使われます。\n", "\n", "```{csv-table}\n", ":header-rows: 1\n", "\n", "\"都市 0\", \"都市 1\", \"都市 2\", \"都市 3\", \"都市 4\"\n", "0, 0, 0, 1, 0\n", "```\n", "\n", "巡回路では都市を 1 回ずつ訪れて最後に最初の都市に戻ってくるので、最初と最後も含めて `NUM_CITIES + 1` 回都市を訪れることになります。したがって、 `(NUM_CITIES + 1) * NUM_CITIES`{l=python} 個のバイナリ変数を用意します。たとえば以下の例は、都市 3 → 都市 1 → 都市 4 → 都市 0 → 都市 2 → 都市 3 の順に 5 つの都市を訪れる巡回路を表現します。\n", "\n", "```{csv-table} バイナリ変数表\n", ":header-rows: 1\n", "\n", "\"都市 0\", \"都市 1\", \"都市 2\", \"都市 3\", \"都市 4\"\n", "0, 0, 0, 1, 0\n", "0, 1, 0, 0, 0\n", "0, 0, 0, 0, 1\n", "1, 0, 0, 0, 0\n", "0, 0, 1, 0, 0\n", "0, 0, 0, 1, 0\n", "```\n", "\n", "(tsp-formulation-constraint)=\n", "\n", "### 制約条件\n", "\n", "上の「バイナリ変数表」において、`(NUM_CITIES + 1) * NUM_CITIES`{l=python} 個のバイナリ変数が好き勝手に 0 または 1 の値をとっても巡回路の表現になるとは限りません。たとえば、すべてのバイナリ変数の値が 1 であるとき、何の巡回路も表しません。したがって、バイナリ変数に制約を課す必要があります。必要な制約は以下の 3 種類です。\n", "\n", "1. バイナリ変数表の最初の行と最後の行は、同じ値をとる。\n", "2. バイナリ変数表の各行のうち、値が 1 となっている箇所はちょうど 1 つだけである。\n", "3. バイナリ変数表の各列において、値が 1 となっている箇所は最後の行を除きちょうど 1 つだけである。\n", "\n", "逆に、以上の 3 つの制約をみたしているバイナリ変数表があるとき、それはある巡回路の表現となっています。\n", "\n", "(tsp-formulation-objective)=\n", "\n", "### 目的関数\n", "\n", "巡回セールスマン問題の目的は、巡回路のうち最も経路長が短いものを探すことです。したがって、目的関数は巡回路の経路長となります。\n", "\n", "都市 i と都市 j の距離が `d[i, j]` であるとします。現在のゴールは、バイナリ変数表の値から巡回路の経路長を求めることです。\n", "\n", "バイナリ変数表の値が分かっている場合、巡回路の経路長は以下の疑似コードのようにして求められます。ただし、`k` 行 `i` 列のバイナリ変数が `q[k, i]` と書けるとします。\n", "\n", "```\n", "route_length = 0\n", "\n", "for k in range(NUM_CITIES):\n", " for i in range(NUM_CITIES):\n", " for j in range(NUM_CITIES):\n", " route_length += (d[i, j] if q[k, i] == 1 and q[k+1, j] == 1 else 0)\n", "```\n", "\n", "しかし Amplify で定式化を行う場合、それぞれの `q[k, i]` はバイナリ変数であることが分かっているのみでまだ何の値を取るかは分からないので、`if q[k, i] == 1 and q[k+1, j] == 1`{l=python} という条件式の判定を行うことはできません。バイナリ変数や数値の足し算および掛け算で目的関数を表す必要があります。\n", "\n", "ここで、2 つのバイナリ変数の積 `q[k, i] * q[k+1, j]` は `q[k, i] == 1 and q[k+1, j] == 1`{l=python} が True のとき 1 となり、False のときは 0 となることを思い出します。すると、目的関数は\n", "\n", "```\n", "route_length = 0\n", "\n", "for k in range(NUM_CITIES):\n", " for i in range(NUM_CITIES):\n", " for j in range(NUM_CITIES):\n", " route_length += d[i, j] * q[k, i] * q[k+1, j]\n", "```\n", "\n", "と書くことができることが分かります。\n", "\n", "### 定式化\n", "\n", "以上の議論を数式で書き下すと、巡回セールスマン問題の定式化は `(NUM_CITIES + 1) * NUM_CITIES` 個のバイナリ変数を用いて\n", "\n", "```{math}\n", "\\begin{align}\n", "&\\text{minimize} \\quad & \\sum_{0 \\leq i, j, k < N} d_{i, j} q_{k, i} q_{k+1, j} & \\\\\n", "&\\text{subject to} \\quad & q_{N - 1, i} = q_{0, i} \\quad & \\text{for} \\quad 0 \\leq i < N, \\\\\n", "& & \\sum_{0 \\leq i} q_{k, i} = 1 \\quad & \\text{for} \\quad 0 \\leq k < N, \\\\\n", "& & \\sum_{0 \\leq k} q_{k, i} = 1 \\quad & \\text{for} \\quad 0 \\leq i < N\n", "\\end{align}\n", "```\n", "\n", "となります。ただし表記の都合上、都市数 `NUM_CITIES` を $N$ と表記しています。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 問題の作成\n", "\n", "まず、問題を作成します。都市数を `NUM_CITIES` として、それぞれの都市の x 座標と y 座標を 0~100 の整数の中からランダムに選びます。以下では、都市の名前を都市 0, 都市 1, ..., 都市 `NUM_CITIES - 1` と呼ぶことにします。\n", "\n", "本サンプルコードでは、簡単のため都市数を 5 とします。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "remove-cell" ] }, "outputs": [], "source": [ "np.set_printoptions(precision=3)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "NUM_CITIES = 5" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "rng = np.random.default_rng()\n", "\n", "x = rng.integers(0, 100, NUM_CITIES)\n", "y = rng.integers(0, 100, NUM_CITIES)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "次に、都市間の距離を表す距離行列を作成します。距離はユークリッド距離を用います。出力は `n * n`{l=python} の NumPy 配列で、`i` 行 `j` 列の要素は都市 `i` と都市 `j` 間の距離、すなわち `(x[i] - x[j]) ** 2 + (y[i] - y[j] ** 2) ** 0.5`{l=python} を表します。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "distance = (\n", " (x[:, np.newaxis] - x[np.newaxis, :]) ** 2\n", " + (y[:, np.newaxis] - y[np.newaxis, :]) ** 2\n", ") ** 0.5\n", "\n", "print(distance)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Amplify SDK による定式化\n", "\n", "Amplify SDK を用いた定式化を行います。\n", "\n", "### 変数の生成 [{octicon}`link`](variables.md)\n", "\n", "まず、決定変数を作成します。`(NUM_CITIES + 1) * NUM_CITIES`{l=python} 個のバイナリ変数が必要です。Amplify SDK では、{py:class}`~amplify.VariableGenerator` を使用することにより 2 次元配列の形で変数を生成することができます。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from amplify import VariableGenerator\n", "\n", "gen = VariableGenerator()\n", "q = gen.array(\"Binary\", (NUM_CITIES + 1, NUM_CITIES))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "バイナリ変数 `q` を表示してみると、以下のようになります。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "q" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(tsp-constraint)=\n", "\n", "### 制約条件の作成 [{octicon}`link`](constraint.md)\n", "\n", "次に、制約条件を作成します。[前述](#tsp-formulation-constraint)の通り、作成するべき制約条件は以下の 3 種類です。\n", "\n", "1. バイナリ変数表の最初の行と最後の行は、同じ値をとる。\n", "2. バイナリ変数表の各行のうち、値が 1 となっている箇所はちょうど 1 つだけである。\n", "3. バイナリ変数表の各列において、値が 1 となっている箇所は最後の行を除きちょうど 1 つだけである。\n", "\n", "#### 1 つ目の制約条件\n", "\n", "1 つ目の制約条件「バイナリ変数表の最初の行と最後の行は、同じ値をとる。」を作成します。最後の行を最初の行と同じ値に固定したいので、`q` の最後の行に最初の行を代入します。`q` の最初の行は `q[0]`、`q` の最後の行は `q[-1]` と書くことができます。\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "q[-1] = q[0]\n", "\n", "q" ] }, { "cell_type": "markdown", "metadata": {}, "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "これで 1 つ目の制約条件を表現することができました。\n", "\n", "```{attention}\n", "このように代入により変数を固定する操作を行う場合、目的関数や他の制約条件を構築する前に行ってください。\n", "```\n", "\n", "これ以降の制約条件の作成においては、バイナリ変数配列 `q` の最後の行のことは考える必要がないため、`q` の上から `NUM_CITIES` 行を切り出した配列を作成しておきます。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "q_upper = q[:-1]\n", "\n", "q_upper" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 2 つ目の制約条件\n", "\n", "2 つ目の制約条件「バイナリ変数表の各行のうち、値が 1 となっている箇所はちょうど 1 つだけである。」を作成します。いくつかのバイナリ変数のうちちょうど 1 つだけが 1 であるという制約条件を課したい場合、{py:func}`~amplify.one_hot` 関数を使います。変数配列の各行に対して one-hot 制約を課したい場合、{py:func}`~amplify.one_hot` 関数の `axis` キーワード引数に 1 を与えます。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from amplify import one_hot\n", "\n", "constraints2 = one_hot(q_upper, axis=1)\n", "\n", "constraints2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`NUM_CITIES` (= 5) 個の制約条件が一気に生成されました。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 3 つ目の制約条件\n", "\n", "3 つ目の制約条件「バイナリ変数表の各列において、値が 1 となっている箇所は最後の行を除きちょうど 1 つだけである。」を作成します。2 つ目の制約条件と同様に、{py:func}`~amplify.one_hot` 関数を使います。変数配列の各列に対して one-hot 制約を課したい場合、{py:func}`~amplify.one_hot` 関数の `axis` キーワード引数に 0 を与えます。\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "constraints3 = one_hot(q_upper, axis=0)\n", "\n", "constraints3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(tsp-objective)=\n", "\n", "### 目的関数の作成 [{octicon}`link`](objective.md)\n", "\n", "次に、目的関数を作成します。[前述](#tsp-formulation-objective)の通り、目的関数は巡回路の経路長であり、以下のように書くことができます。\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "route_length = 0\n", "\n", "for k in range(NUM_CITIES):\n", " for i in range(NUM_CITIES):\n", " for j in range(NUM_CITIES):\n", " route_length += distance[i, j] * q[k, i] * q[k + 1, j]\n", "\n", "print(route_length)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "あるいは、以下のように書いても同じ効果があります。こちらの方が Amplify SDK の機能をフルに活用しているため高速で、NumPy ライブラリに馴染みがある方には特におすすめです。詳細は[](optimization.md)を参照してください。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from amplify import Poly, einsum\n", "\n", "q1 = q[:-1]\n", "q2 = q[1:]\n", "route_length: Poly = einsum(\"ij,ki,kj->\", distance, q1, q2) # type: ignore\n", "\n", "print(route_length)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(tsp-model)=\n", "\n", "### 組合せ最適化モデルの作成 [{octicon}`link`](model.md)\n", "\n", "[制約条件の作成](#tsp-constraint)、[目的関数の作成](#tsp-objective) で作成した目的関数と制約条件を組み合わせて、組合せ最適化のモデルを作成します。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model = route_length + (constraints2 + constraints3) * np.max(distance)\n", "\n", "print(model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "制約条件に距離行列の最大値をかけているのは、制約条件に重みをつけるためです。今回使用するソルバーである Amplify AE においては、制約条件に適切な重みを指定しないと、ソルバーが制約条件をみたそうとするよりも目的関数を小さくする方向に動いてしまい、良い解を発見することができなくなってしまいます。詳しくは [](penalty.md) を参照してください。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## ソルバークライアントの作成 [{octicon}`link`](clients.md)\n", "\n", "ソルバークライアントを作成し、求解を行うためのソルバーの指定およびソルバーのパラメータの設定を行います。Amplify SDK はさまざまなソルバーに対応していますが、本サンプルコードでは Amplify AE を使用します。Amplify AE に対応するソルバークライアントクラスは {py:class}`~amplify.FixstarsClient` クラスです。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from amplify import FixstarsClient\n", "\n", "client = FixstarsClient()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Amplify AE の実行に必要な API トークンを設定します。\n", "\n", "```{tip}\n", "[ユーザ登録](https://amplify.fixstars.com/register)を行うと、評価・検証目的に使える API トークンを無料で入手できます。\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "skip-execution" ] }, "outputs": [], "source": [ "client.token = \"xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx\"" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "remove-cell" ] }, "outputs": [], "source": [ "client.token = \"\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "ソルバーのタイムアウト値を設定します。Amplify AE のタイムアウト値の単位は ms ですが、タイムアウト値の単位はソルバーによって様々なため、{py:mod}`datetime` モジュールを使用して統一的にタイムアウト値を指定できます。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import datetime\n", "\n", "client.parameters.timeout = datetime.timedelta(seconds=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "これでソルバーの設定は完了です。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 求解の実行 [{octicon}`link`](solve.md)\n", "\n", "作成した組合せ最適化モデルとソルバークライアントを使用してソルバーの実行を行い、巡回セールスマン問題の解を求めます。" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "from amplify import solve\n", "\n", "result = solve(model, client)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "目的関数の値、つまり巡回路の長さは以下のように表示できます。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "result.best.objective" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "最適解における変数の値は、NumPy の多次元配列の形式で以下のようにして取得できます。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "q_values = q.evaluate(result.best.values)\n", "\n", "print(q_values)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 結果の確認\n", "\n", "Amplify SDK のチュートリアルとしては以上で十分ですが、結果が正しそうであることを可視化しておきます。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "まず、都市の位置を可視化します。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.scatter(x, y)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "次に、巡回路を可視化します。`q_values` は置換行列 (の最後に 1 行足したもの) と見なせるので、各都市を順に並べたベクトルに `q_values` を左から掛けることにより、巡回路で訪問する順に並び替えることができます。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "route_x = q_values @ x # 移動ルートの x 座標を取得\n", "route_y = q_values @ y # 移動ルートの y 座標を取得\n", "\n", "plt.scatter(x, y)\n", "plt.plot(route_x, route_y)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "最短路が求められていることが分かります。" ] } ], "metadata": { "kernelspec": { "display_name": "amplify-lrEeNlW4", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.12" } }, "nbformat": 4, "nbformat_minor": 2 }