「Pythonを用いた巡回セールスマン問題の実装、解決、および可視化」

「Pythonを使った巡回セールスマン問題の実装、解決、そして可視化の方法」

数学からPythonへの最適化モデルの翻訳、最適化、および解の可視化によるモデリングエラーの迅速なフィードバックの獲得

👁️ これは、“Pythonにおける観光のための知的な意思決定支援システム”プロジェクトをカバーするシリーズの第3回記事です。 全体的なプロジェクトの概要を把握するために、ぜひチェックしてみてください。もしTSPのモデルをPythonで実装する方法にしか興味がなければ、この記事も適しています。自己完結型な記事であり、依存関係のインストール、分析、コード、結果の解釈、およびモデルのデバッグのすべての手順を詳しく解説しています。

この記事では、スプリント2で停止した箇所から旅を続けます。ここでは、前の記事で定式化した数学的なモデルをPythonで実装し、Pyomoを使用してベストプラクティスに従います。その後、モデルを最適化し、その解を視覚化し分析します。教育のために、初期のモデルの定式化が不完全であることがわかったため、定式化を修正するために最初の原理から派生できる制約条件を示します。これらの新しい制約条件はPyomoモデルに追加され、新しい解が再度分析され検証されます。

目次

1. Pyomoを使用してPythonでモデルを実装する

• 1.1. 依存関係のインストール
• 1.2. 数学からコードへ

2. モデルの解決と検証

• 2.1. モデルの解決
• 2.2. 結果の視覚化
• 2.3. 結果の分析

3. 定式化の修正

• 3.1. 動機づけのアイデア
• 3.2. 論理的な含意として動機づけのアイデアを表現する
• 3.3. 線形制約として論理的な含意を定式化する
• 3.4. “big M”の適切な値を推論する

4. 新しい定式化の実装と検証

• 4.1. Pyomoモデルの拡張
• 4.2. 更新されたモデルの解のプロット

5. 結論（次のスプリントについて）

📖 前の記事を読んでいなくても大丈夫です。数学的な定式化もここで明記されています（ただし、導出はされていません）。各モデルコンポーネントは、それに対応するコード実装と一緒に示されています。もし何がどこから来て何を意味するのか理解できない場合は、“スプリント2”の記事を読んでみてください。さらなる背景情報が必要な場合は、“スプリント1”の記事も読んでみてください。

1. Pyomoを使用してPythonでモデルを実装する

データサイエンスで最も優れた言語であるため、Pythonが使用されており、大規模なモデルに効果的に対処するための最高の（オープンソースの）ライブラリであるPyomoが使用されています。

• このAI研究では、優れた画像分類のためにランダムスライスデータ拡張（RSMDA）を提案します：ニューラルネットワークの精度と堅牢性を向上させるための新しいアプローチ
• 「CNNによる特徴抽出の探求」
• 「それに関する長いものと短いもの：ドキュメントの意味論を端から端まで捉えるための比例関係に基づく関連性」

このセクションでは、定式化で定義された各モデルコンポーネントを逐一説明し、それがPyomoコードにどのように変換されるかを説明します。間違いがないように心がけましたが、違う思いがある場合は、コメントで質問を残してください。

免責事項：対象読者はPyomoやモデリングに新しいことを想定しており、認知負荷を軽減するために、簡潔でシンプルな実装がプログラミングのベストプラクティスよりも優先されます。今の目標は最適化モデリングを教えることであり、ソフトウェアエンジニアリングではありません。この概念実証がより洗練されたMVPに進化するにつれて、コードは将来のイテレーションで段階的に改善されます。

1.1. 依存関係のインストール

時間がない人へ

pyomo、networkx、pandasのライブラリとglpkパッケージをインストール（または既にインストールされていること）するようにしてください。

📝 glpkパッケージにはGLPKソルバーが含まれており、モデリングした問題を最適化するために使用します。Pyomoを使用して問題のモデルを作成し、それをGLPKに渡して最適化プロセス自体を実行するアルゴリズムを実行します。その後、GLPKはPyomoモデルオブジェクトに対して解を返し、それらはモデルの属性として格納されるため、Pythonを離れることなく便利に使用できます。

GLPKをインストールする推奨方法は、PyomoがGLPKソルバーを簡単に見つけられるようにするために、condaを介して行うことです。一括ですべての依存関係をインストールするには、次のコマンドを実行してください：

conda install -y -c conda-forge pyomo=6.5.0 pandas networkx glpk

整理された人へ

このシリーズの記事に従うために必要なすべてのライブラリをインストールするために、独立した仮想環境を作成することをお勧めします。以下のテキストをコピーして

name: ttp  # 旅行案内人問題channels:  - conda-forgedependencies:  - python=3.9  - pyomo=6.5.0  - pandas  - networkx  - glpk  # 外部ソルバー（最適化モデルに使用）  - jupyterlab  # IDEとしてJupyter Labを使用しない場合は、この行をコメントアウト

そして、それをenvironment.ymlという名前のYAMLファイルに保存してください。同じ場所でAnacondaプロンプトを開き、次のコマンドを実行します：

conda env create --file environment.yml

数分後、環境が作成され、すべての依存関係が内部にインストールされます。環境に「入る」には、conda activate ttpを実行し、ターミナルでjupyter labを実行してJupyter Labを起動し、コーディングを始めます。

1.2. 数式がコードに変わる

まず、GLPKソルバーがPyomoに見つけられるようになっているか確認しましょう

### =====  コードブロック 3.1  ===== ###import pandas as pdimport pyomo.environ as pyoimport pyomo.versionimport networkx as nxsolver = pyo.SolverFactory("glpk")solver_available = solver.available(exception_flag=False)print(f"Solver '{solver.name}' available: {solver_available}")if solver_available:    print(f"Solver version: {solver.version()}")print("pyomo version:", pyomo.version.__version__)print("networkx version:", nx.__version__)

Solver 'glpk' available: TrueSolver version: (5, 0, 0, 0)pyomo version: 6.5networkx version: 2.8.4

⛔ メッセージ'glpk' available: Falseが表示された場合、ソルバーが正しくインストールされませんでした。以下のいずれかを試して問題を修正してください：

– インストール手順を注意深く再実行する

– ベース（デフォルト）環境でconda install -y -c conda-forge glpkを実行する

– 動作する別のソルバーをインストールしてみる

その後、距離データのCSVファイルを読み取ります

### =====  コードブロック3.2  ===== ###path_data = (    "https://raw.githubusercontent.com/carlosjuribe/"    "traveling-tourist-problem/main/"    "Paris_sites_spherical_distance_matrix.csv")df_distances = pd.read_csv(path_data, index_col="site")df_distances

そして、[敏捷なオペレーション研究ワークフロー]の「ステージ4」に入ります。以下の緑色のブロックで示されています:

タスクは、以前に作成した数学モデルをそのままコードに実装することです。

👁️ モデルの実装を容易にするために、Pythonオブジェクトを作成することは許可されていますが、コーディング中に基礎となるモデルを修正することは許可されていません。これにより、数学モデルとコンピュータモデルが同期しなくなるため、後でモデルのデバッグが非常に困難になります。

空のPyomoモデルを作成し、モデルのコンポーネントを属性として保存します:

model = pyo.ConcreteModel("TSP")

1.2.1. セット

サイトのセット𝕊 = {Louvre、エッフェル塔、…、ホテル}を作成するために、データフレームのインデックスからサイトの名前を取得し、Pyomo Setであるsitesとして使用します:

### =====  コードブロック3.3  ===== ###list_of_sites = df_distances.index.tolist()model.sites = pyo.Set(initialize=list_of_sites,                       domain=pyo.Any,                       doc="訪れるすべてのサイトのセット（𝕊）")

派生セットの作成

フィルター𝑖 ≠ 𝑗を構築ルール（Python関数_rule_domain_arcs）に格納し、Setの初期化時にこのルールをfilterキーワードに渡します。このフィルターはサイト（𝕊 × 𝕊）の直積に適用され、ルールを満たさない直積のメンバーはフィルタリングされます。

### =====  コードブロック3.4  ===== ###def _rule_domain_arcs(model, i, j):    """ サイトを接続する可能なすべてのアーク """    # サイトiとサイトjが異なる場合のみペア（i、j）を作成する    return (i, j) if i != j else None rule = _rule_domain_arcsmodel.valid_arcs = pyo.Set(    initialize=model.sites * model.sites,  # 𝕊 × 𝕊    filter=rule, doc=rule.__doc__)

1.2.2. パラメータ

パラメータ

𝐷ᵢⱼ ≔ サイト𝑖とサイト𝑗の距離

はコンストラクタpyo.Paramを使用して作成されます。先頭の（位置指定の）引数はドメイン𝔸(model.valid_arcs)であり、キーワード引数initializeには別の構築ルール_rule_distance_between_sitesがあります。このルールは、𝑖、𝑗∈𝔸の各ペアに対して評価され、評価ごとに距離の数値がデータフレームdf_distancesから取得され、内部的にペア（𝑖、𝑗）とリンクされます:

### =====  コードブロック3.5  ===== ###def _rule_distance_between_sites(model, i, j):    """ サイトiとサイトjの距離（𝐷𝑖𝑗） """    return df_distances.at[i, j]  # データフレームから距離を取得するrule = _rule_distance_between_sitesmodel.distance_ij = pyo.Param(model.valid_arcs,                               initialize=rule,                               doc=rule.__doc__)

1.2.3. 決定変数

𝛿ᵢⱼが𝐷ᵢⱼと同じ「インデックスドメイン」を持っているため、このコンポーネントを構築する方法は非常に似ていますが、ここでは建設ルールは必要ありません。

model.delta_ij = pyo.Var(model.valid_arcs, within=pyo.Binary,                          doc="サイトiからサイトjへ移動するかどうか（𝛿𝑖𝑗）")

pyo.Varの最初の位置引数は、そのインデックスセット𝔸に予約されており、変数の「型」はキーワード引数withinで指定されます。ここで「変数の型」とは、変数が取り得る値の範囲のことを指します。ここでは、𝛿ᵢⱼの範囲は0と1のみなので、バイナリの型です。数学的には、𝛿ᵢⱼ ∈ {0, 1}と書くことになるのですが、これを示すために別々の制約を作成するのではなく、直接Pyomoで指定することができます。within=pyo.Binaryを変数の作成時に設定することで、直接示すことができます。

1.2.4. 目的関数

目的関数を構築するためには、式を関数に「格納」することができます。この関数は、式を構築するために必要なモデルコンポーネントを取得するために使用されるモデル引数のみを受け取ります。

### =====  コードブロック3.6  ===== ###def _rule_total_distance_traveled(model):    """ 総移動距離 """    return pyo.summation(model.distance_ij, model.delta_ij)rule = _rule_total_distance_traveledmodel.obj_total_distance = pyo.Objective(rule=rule,                                          sense=pyo.minimize,                                          doc=rule.__doc__)

数学的な式と関数のreturn文との類似性に注目してください。最小化問題であることをsenseキーワードで指定しています。

1.2.5. 制約条件

前の記事で説明したように、各サイトが一度だけ訪れられるようにする便利な方法は、各サイトが一度だけ「入る」と同時に一度だけ「出る」ように強制することです。

各サイトは一度だけ入る

def _rule_site_is_entered_once(model, j):    """ 各サイトjは他のサイトからちょうど1回訪れられる必要がある """    return sum(model.delta_ij[i, j] for i in model.sites if i != j) == 1rule = _rule_site_is_entered_oncemodel.constr_each_site_is_entered_once = pyo.Constraint(                                          model.sites,                                           rule=rule,                                           doc=rule.__doc__)

各サイトは一度だけ退出します

def _rule_site_is_exited_once(model, i):
    """ サイト i は他のサイトに正確に一度退出する必要があります """
    return sum(model.delta_ij[i, j] for j in model.sites if j != i) == 1

rule = _rule_site_is_exited_oncemodel.constr_each_site_is_exited_once = pyo.Constraint(
                                          model.sites,
                                          rule=rule,
                                          doc=rule.__doc__
)

1.2.6. モデルの最終検査

モデルの実装は完了しました。全体を確認するためには、model.pprint() を実行し、宣言の見落としや誤りがないかを少し探索します。

モデルのサイズを把握するために、変数の数と制約の数を表示します：

def print_model_info(model):
    print(f"Name: {model.name}",
          f"変数の数: {model.nvariables()}",
          f"制約の数: {model.nconstraints()}",
          sep="\n- ")

print_model_info(model)#[Out]:# Name: TSP
#  - 変数の数: 72
#  - 制約の数: 18

制約や変数の数が100以下であれば、この問題は小規模であり、ソルバーによって比較的高速に最適化されます。

2. モデルの解決と検証

2.1. モデルの解決

[AORフローチャート]の次のステップは、モデルを最適化し、解を検証することです：

res = solver.solve(model)  # モデルを最適化する
print(f"最適解が見つかりました: {pyo.check_optimal_termination(res)}")
# [Out]: 最適解が見つかりました: True

よかったですね！ソルバーがこの問題の最適解を見つけました！パリに到着したときにどのツアーを辿るかを知るために、それをチェックしましょう！

非常に簡単なチェックのために、model.delta_ij.pprint() を実行することができます。これは、𝛿ᵢⱼ変数の（最適な）値すべてを、0または1であるかのいずれかとして表示します：

選択されたアークのリストを見てツアーを直感的に視覚化することは困難であり、さらにはモデルの正確な定式化を検証するために分析することも困難です。

解を本当に理解するためには、それを視覚化する必要があります。

2.2. 結果の視覚化

画像は千の価値がある

ノードとアークを扱っているため、解をグラフとしてプロットするのが最も簡単な方法です。これは概念実証であり、素早く効果的なフィードバックが美しさに優るものです。洞察力のある視覚化は、動作するMVPがあるまで待つことができます。今のところ、解を効率的にプロットするためのヘルパー関数をいくつか作成しましょう。その

extract_solution_as_arcs 関数は解かれたPyomoモデルを受け取り、解から「選択されたアーク」を抽出します。次に、plot_arcs_as_graph 関数はアクティブなアークのリストをグラフオブジェクトに保存して、分析を容易にし、ホテルが唯一の赤いノードであるようにそのグラフをプロットします。最後に、plot_solution_as_graph 関数は上記の二つの関数を呼び出して、与えられたモデルの解をグラフとして表示します。

### =====  コードブロック 3.7  ===== ###def extract_solution_as_arcs(model):    """ 解かれたモデルからアクティブ（選択された）アークのリストを抽出するため、    1のバイナリ変数を持つアークのみを残す """    active_arcs = [(i, j)                   for i, j in model.valid_arcs                   if model.delta_ij[i, j].value == 1]    return active_arcsdef plot_arcs_as_graph(tour_as_arcs):    """ アークのリストを受け取り、networkxグラフに変換し、描画する """    G = nx.DiGraph()    G.add_edges_from(tour_as_arcs)  # グラフとしてソリューションを保存    node_colors = ['red' if node == 'hotel' else 'skyblue'                    for node in G.nodes()]    nx.draw(G, node_color=node_colors, with_labels=True, font_size=6,             node_shape='o', arrowsize=5, style='solid')    def plot_solution_as_graph(model):    """ モデルの解をグラフとして描画する """    print(f"総距離: {model.obj_total_distance()}")        active_arcs = extract_solution_as_arcs(model)    plot_arcs_as_graph(active_arcs)

さて、それでは解決策がどのように見えるかを見てみましょう：

plot_solution_as_graph(model)

さて、これは明らかに私たちの期待とは異なる結果です！全ての場所を訪れるはいいですが、ただ小さな切れ切れのスタート間所が訪れられるだけです。技術的には、指定した制約はきちんと守られています：各場所は一度だけ入場し、一度だけ出場するが、全体的な結果は私たちが意図したような一つのツアーではなく、いくつかのサブツアーのグループです。つまり、前の記事で行った仮定は間違っており、解が許可されないことをエンコードするモデルに不足しているものがあることを意味します。

2.3. 結果の分析

何が間違ってしまったのでしょうか？

モデルの解が意味をなさない場合、唯一の可能な説明は、モデルが間違っていることです¹。

ソルバーは、私たちが解きたい問題にマップしないモデルを与えましたが、それは確かにモデルの最適解です。課題は、なぜ間違ったのか、どこでミスをしたのかを見つけ出すことです。反省して考えると、「4.4. 制約の作成」セクションの最後の2段落で私たちが行った疑わしい仮定を指摘するのは明らかです。数学モデルを設計する際に、私たちは（現在は誤りであることを知っている）「単一のツアーの形成は、各場所が一度だけ訪れるという2つの制約から自然に続く」と仮定しました。しかし、私たちが視覚化したように、それは起こりません。なぜでしょうか？

エラーの根本原因は、「口に出さない常識」にあります：モデルに明示することを忘れてしまうほど明らかな世界に関する知識

私たちは、サイトを訪れる際にテレポーテーションが不可能であることを暗黙のうちに知っていましたが、モデルに明示的に伝えていませんでした。そのため、いくつかのサイトを接続している小さなサブツアーが観察され、すべてのサイトが繋がっていないのです。「モデルは」、サイトから別のサイトにテレポートすることは許可されていると「思っています」、ただし、サイトに着いたら一度だけ出発し、一度だけ到着するだけです（図3.2を再確認してください）。サブツアーを見るのは、モデルにツアーの距離を最小化するよう指示したからであり、たまたまテレポートが距離を節約するために役立つからです。

<!–したがって、これらの複数の小旅行の形成を防止する必要があります 、現実的な解決策を得るためには。 いくつかの新しい制約を設計する必要があります 、その制約をモデルに伝えるもしくは同等の解決策は単一のツアーでなければならない ことを 示す方法があります。 後者を選びましょう。また、 初歩の原理から導き出す 直感的でうまく機能する制約セットを1つのセット作りましょう。

3. フォーミュレーションの修正

[アジャイルオペレーションリサーチのワークフロー]に従って、私たちは現在モデルの再構成フェーズにいます。モデルの再構成は、改善するか修正するかについてのものです。私たちの場合は修正についてです。

私たちが望んでいることはわかっています。解決策を初めてのサイトで始まり、ホテルで終わる単一のツアーにすることです。 課題は、その要件を一連の線形制約にエンコードする方法です 。以下は、ツアーの性質から導かれるアイデアの一つです。

3.1. 動機づけのアイデア

ホテルを含む𝑁個のサイトを巡る必要があります。ホテルで始まるため、訪れるサイトは𝑁−1個です。もし、それらのサイトが「時系列の順序」で訪れた順序を追跡し、最初の目的地（ホテルの後）には番号1が与えられ、2番目の目的地には番号2が与えられる、という風に続いていくとすると、ホテルに戻る前の最後の目的地は番号𝑁−1が与えられます。この訪問順序を追跡するために使われるこれらの番号を「ランク」と呼びます。したがって、ツアー内で前のサイトに比べると、（ホテル以外の）サイトのランクは常に1高いことになります。もし、任意の解がこのようなパターンを実現するような制約セットを提示できたら、モデルに「時系列の順序」の要件を緩く導入することになります。そして、実際にそうできることがわかります。

3.2. 動機づけのアイデアを論理的な含意として表現する

💡任意の解に満たすことを望む「パターン」です:

ホテルを除く任意のサイトのランクは、その前のサイトのランクよりも必ず1高くなければなりません

このパターンを 論理的な含意 として再表現することができます。次のように：「ホテル 𝐻 を含まないすべての弧(𝑖,𝑗)について、もし𝑗が𝑖の直後に訪れた場合は、サイト𝑖のランクよりもサイト𝑗のランクは必ず1高くなければなりません」と表現されます。この記述は数学的に次のように表されます：

ここで、𝑟ᵢは訪問順序の（まだ未知の）追跡に必要な新しい変数です。これらを決定変数と区別するために、「ランク変数」と呼びます。右辺は「すべての𝑖および𝑗がホテルを含まないすべてのサイトの集合に属する」と読みます。表記の便宜のために、新しい集合𝕊* を定義して、ホテルを除くすべてのサイトを格納します（𝐻で表されます）：

これにより、ランク変数を簡潔に定義することができます：

👁️ ホテルには関連するランク変数が存在しないことが重要です。なぜなら、ホテルはツアーの起点および最終目的地であり、これは「ツアー内のランク変数が常に増加する」というパターンを破る条件です。したがって、各サイトのランクは常に新しいアークが取られるごとに増加することが強制され、閉ループは唯一のランク変数を持たない唯一のサイトで閉じる場合を除き禁止されます：ホテル

𝑟ᵢの上限はその説明から導かれます：ランクは1から始まり、𝕊 *内のすべてのサイトが訪れられるまで単調に増加し、したがって | 𝕊 * |（非ホテルサイトのセットのサイズ）で終了します。また、正の実数の任意の値を取ることができるようにします：

3.3. 論理的含意を線形制約として表現する

今の課題は、この論理的含意を線形制約のセットに変換することです。幸いなことに、線形不等式は論理的な含意を強制するためにも使用されます。有限なリソースの制約だけでなく。

これを行う一つの方法は、その制約を宣言するためのビッグM法と呼ばれる方法で、この方法では、気になる条件が満たされた場合、制約が適用され（有効になり）、気になる条件が満たされない場合、制約が冗長になり（無効になる）ように制約を宣言します。このテクニックは「ビッグM」と呼ばれます。なぜなら、制約に現れる定数値𝑀は、各ケースで無効にするために十分に大きくなければならないからです。制約に𝑀が現れないときは、制約が「有効」であり、望ましい含意を強制しています。

しかし、制約が「有効」かどうかは何によって決まるのでしょうか？ 短い答えは、制約が適用される決定変数の値そのものです。それがどのように機能するか見てみましょう。

望ましい含意は、𝑟ⱼ = 𝑟ᵢ + 1 が 𝛿ᵢⱼ = 1 の場合にのみ成立することです。この式において1を𝛿ᵢⱼで置き換えることができ、これにより

これは、𝛿ᵢⱼ = 1 の場合に成立してほしい関係ですが、𝛿ᵢⱼ = 0 の場合には成立してほしくありません。 𝛿ᵢⱼ = 0 の場合を「修正」するために、冗長項𝑅𝑇を追加します。この冗長項の役割は、𝛿ᵢⱼ = 0 の場合にのみ制約を「非活性化すること」です。したがって、この冗長項には変数𝛿ᵢⱼを含める必要があります。

この文脈では、「制約の非活性化」とは、「制約が冗長になる」という意味です。なぜなら、冗長な制約はモデル内の存在しない制約と同じ効果を持つからです。

RTの式を導く方法を見てみましょう。RTの式は、以下の特性を満たす必要があります：

点(1)を満たすためには、𝑅𝑇(𝛿ᵢⱼ = 1) = 0が必要です。したがって、𝑅𝑇の式は乗数(1 – 𝛿ᵢⱼ)を含まなければなりません。なぜなら、𝛿ᵢⱼが1の場合には0になるからです。この形式では、𝑅𝑇は𝛿ᵢⱼが1の場合に「消える」または𝛿ᵢⱼが0の場合には定数（Mと呼びましょう）に「縮小される」ようになります。したがって、冗長性項の候補は次のようになります。

ここで、𝑀は問題データから決定する必要があります（詳細は後で述べます）。

（2）をすべての可能な𝑖と𝑗に対して満たすためには、式（3.11）の等号を不等号（=を≥に）にする必要があり、定数𝑀を見つける必要があります。この定数𝑀は、𝑟ⱼと𝑟ᵢがどのような値を取るにしても常に制約を満たすために絶対値が十分に大きいものである必要があります。これが「big M」の「big」の由来です。

十分に大きな定数𝑀を見つけた後、私たちの「論理的な含意」制約は次のようになります。

これらの制約をモデルに導入することで、解が単一のツアーになるようにできるでしょう。ただし、制約が望んだ効果を持つためには、Mの適切な値を最初に指定する必要があります。

3.4. 「big M」の適切な値を導く

𝑅𝑇(𝛿ᵢⱼ = 0) = 𝑀とすることを目標としているため、式（3.12）の一般的な制約に𝛿ᵢⱼ = 0を代入することによって、𝑀の適切な値を導くことができます。

すべての非ホテルサイトの𝑖、𝑗に対して𝑟ⱼ − 𝑟ᵢ ≥ 𝑀が満たされるために、𝑟ⱼ − 𝑟ᵢの下限が𝑀よりも大きくなる必要があります。𝑟ⱼ − 𝑟ᵢの下限（LB）は、𝑟ⱼ − 𝑟ᵢが取り得る最小値であり、次の式によって得ることができます。

𝑟ᵐⁱⁿが最小のランクであり、𝑟ᵐᵃˣが最大のランクである場合、不等式（1）が全てのサイトの全てのランクに対して真であるためには、以下の不等式も成り立たなければなりません:

この不等式により、big-M法が機能するために𝑀が取る必要がある最小の値がわかります。

𝑟ᵐⁱⁿと𝑟ᵐᵃˣの値は何ですか？ 慣例では、最初の訪問サイトにランク1を与えています（つまり、𝑟ᵐⁱⁿ = 1）。ランクはそれぞれの訪問サイトで1ユニットずつ増加するため、ツアーで最後のホテル以外のサイトは最大のランクを持ちます。非ホテルのサイトの数は変化するため、一般的な表現が必要です。 𝕊*セットはすべてのホテル以外のサイトを含んでいることを覚えており、求める数は𝕊*セットの要素数です（数学表記では | 𝕊* |）。したがって、𝑟ᵐⁱⁿ = 1 and 𝑟ᵐᵃˣ = | 𝕊* |となります。これを式（3.14）に代入すると、Mの正しい値が得られます。

𝕊*は必ず2つ以上のサイトを訪れる必要があります（そうでなければ、最初から意思決定問題は存在しない）。したがって、このモデルでは「big M」の値は常に「マイナスの大きな値」です。

理論的な値と計算上の値

理論的には、この計算を避けるために、ここで導かれたよりも「もっとマイナス」の値を𝑀として任意に選ぶことが許されていますが、これは良いやり方ではありません。 𝑀があまりにも大きくなると（絶対値で）、ソルバのアルゴリズムにパフォーマンスの問題を引き起こしたり、最悪の場合にはソルバが「不可能な解を可能な解と見なすことさえあります。そのため、推奨される方法は、問題のデータから推測されたを得るために、適切かつ十分に大きな値を導出することです。

これで「big M」の適切な値が導かれたので、より簡単に再利用できるように新しいモデルパラメータに保存します。これにより、副巡回制約セットが準備され、その「完全な形式」は次のとおりです:

透視を保つために、これは実際には前述の論理含意の “制約等価” です。

おめでとうございます！私たちはついに、モデルに追加できる制約のセットを持っています。それらを考案するのは難しい部分でした。さあ、モデルに追加することで、実際に巡回セールスマン問題 (subtours) がなくなることを確認しましょう。

4. 新しい公式の実装と検証

4.1. Pyomoモデルへの新しい公式の追加

モデルの見直しと修正では、いくつかのセット、パラメータ、変数、および制約の追加が必要でした。これらの新しいモデルコンポーネントをPyomoモデルに追加していきましょう。初期の定式化フェーズと同じ順序で追加します。

4.1.1. セットおよびパラメータ

• 興味のある場所、𝕊*：ホテルを除くすべての場所の集合：

Pyomo SetオブジェクトはPythonのセットと互換性のある操作を持っているため、直接PyomoセットとPythonセットの差をとることができます：

model.sites_except_hotel = pyo.Set(    initialize=model.sites - {'ホテル'},     domain=model.sites,    doc="興味のある場所、つまりホテル以外のすべての場所（𝕊*）")

• big M、サブツアー除去制約のための新しいパラメータ：

model.M = pyo.Param(initialize=1 - len(model.sites_except_hotel),                    doc="ある制約を冗長にするためのbig M")

4.1.2. 補助変数

• 順位変数、r_i：訪れる場所の順序を追跡するために使用：

model.rank_i = pyo.Var(    model.sites_except_hotel,  # i ∈ 𝕊* (index)    within=pyo.NonNegativeReals,  # rᵢ ∈ ℝ₊ (domain)    bounds=(1, len(model.sites_except_hotel)),  # 1 ≤ rᵢ ≤ |𝕊*|    doc="各場所の順位、訪れる順序を追跡")

コメントでは、変数の完全な数学的定義の要素がPyomo変数宣言関数 pyo.Var の引数にどのようにマップされるかがわかります。Pyomoモデルを構築する前に、明確に定義された数学モデルがあることの価値を理解していただければ幸いです。実装は自然に進行し、エラーが発生する可能性が少なくなります。

4.1.3. 制約

• 解はホテルから出発し、ホテルで終了する単一のツアーでなければなりません：

def _rule_path_is_single_tour(model, i, j):    """ それぞれの非ホテルのサイトの組（i、j）について、
    サイトjがサイトiから訪れられる場合、jのランクはiのランクよりも厳密に大きくなければなりません。"""    if i == j:  # サイトが重なる場合、制約は作成されずスキップします
        return pyo.Constraint.Skip
        
    r_i = model.rank_i[i]
    r_j = model.rank_i[j]
    delta_ij = model.delta_ij[i, j]
    
    return r_j >= r_i + delta_ij + (1 - delta_ij) * model.M# ホテル以外のサイトの直積、制約をインデックス化するための変数
non_hotel_site_pairs = model.sites_except_hotel * model.sites_except_hotel
rule = _rule_path_is_single_tour
model.constr_path_is_single_tour = pyo.Constraint(
    non_hotel_site_pairs,
    rule=rule,
    doc=rule.__doc__)

Pyomoモデルは更新されました。サブツアー除外制約の追加後、どれくらい成長しましたか？

print_model_info(model)# [Out]:# 名前： TSP#   - 変数の数： 80#   - 制約の数： 74

変数の数は72から80に増え、制約の数は18から74に増えました。明らかに、この表現は変数よりも制約の数が多くなるため、通常、モデルをより「現実的」にするための「代価」になります。現実性では、データが変更されない限り、許容解の数を制限することになります。

いつものように、model.pprint()でモデルの構造を確認できます。モデル要素の数が増えると、この「プリント」はその価値を急速に失いますが、モデルが何でできていて、どれほど大きいかについての簡単な概要を提供してくれます。

4.2. 更新されたモデルの解をプロットする

更新されたモデルを解決し、新しい解をプロットしましょう。祈っています。

res = solver.solve(model)  # モデルを最適化する
solution_found = pyo.check_optimal_termination(res)
print(f"最適解が見つかりました：{solution_found}")# [Out]: 最適解が見つかりました：True
if solution_found:
    plot_solution_as_graph(model)

これになりました！これは、新しい制約の追加後に期待していたものです：サブツアーは形成されておらず、ソリューションパスは単一のツアーとなっています。

明らかに、この解決済みモデルの目的関数値は、不完全なモデルで得た5.7 kmではなく、現在は14.9 kmです。

👁️ グラフの描画は地図上のツアーではありません

この画像は、単なるグラフの一つの可能な描画であり、地理的な軌跡ではありません。見ている円やリンクは、地理的な空間上の実際の経路に対応していないことに注意してください（どのようにして実現できますか、なぜなら、それを作成する際に地理情報を使用していないからです）。複数回plot_solution_as_graph(model)を実行すると、ノードは異なる位置を取ります。グラフは、任意の種類のエンティティ間の任意の種類の関係を表す「ポイント」を「リンク」で接続する抽象的な数学的構造です。私たちはここでグラフを使用して、解の妥当性を研究するために使用し、パリの実際のツアーを視覚化するためのものではありません。[この記事]でそれをします。

5. 結論（または次のスプリントの計画）

この解の最終的な検証により、この更新されたバージョンのモデルは旅行セールスマン問題の任意のインスタンスを解決できると結論付けますので、これは成功した概念実証（POC）と見なすことができます。

💡 スプリントごとの解の進化

このPOCは、元の複雑な観光問題を解決するわけではありませんが、より高度な解決策に向けた最初のステップとして提案された最小限の価値のある問題を解決します。したがって、これは証拠に基づく方法で、より複雑な問題のための有益なソリューションに対して証明されている形となります。手元に最小限の動作する例があることで、進むべき方向をより良く評価し、より成熟したバージョンのソリューションが達成されるまで一時的に簡略化できる内容を把握することができます。常に何か有用なものを持つことで、価値のあるシステムを開発し、satisficed になるまで進化する事ができます。それが効果的なソリューションへの「アジャイルの道」の本質です。

このアプローチの妥当性が証明されたことにより、私たちはそれを拡大し、洗練させる必要があります。徐々に問題のより多くの特徴を包含できるようにし、各イテレーションで徐々に価値のある解決策を提供します。このPOCでは、基本モデルの設計と公式に焦点を当てたため、一定のサイトセットとその距離行列を所与とする必要がありました。もちろん、これには制限があり、次のステップは、任意の数のサイトを受け入れるモデルを持つことです。そのためには、サイトのリストと地理座標が与えられた場合に自動的に距離行列を生成する方法が必要です。次のスプリントの目標です。

5.1. 次にやること

次の記事（スプリント4）では、任意のサイトのリストから自動的に距離行列を生成するクラスを作成します。この機能は、ここで構築したモデルと組み合わせると、異なる入力に対して多くのモデルを素早く解決し、比較することができます。また、問題を一般化することは、将来のスプリントでの感度分析とシナリオ分析が容易になります。さらに、この概念実証を「MVPステータス」にアップグレードするにつれて、事物を整理し、拡張性を備えたオブジェクト指向のコードを使用し始めます。流れを無駄にせず、今すぐ進んでみてください！

脚注

1. 実際には、誤った結果の別の原因もあります：モデルがフィードされるデータも誤っている場合です。しかし、ある意味では、モデルを「モデルのインスタンス」すなわち具体的なデータを持つ具体的なモデルと見なすと、データが間違っている場合はもちろんモデルも間違っています。これが私の意図したものです。

読んでいただき、次回をお楽しみに！📈😊

気軽にフォローしていただいたり、質問したり、コメントでフィードバックをくれたり、LinkedInで連絡してください。

We will continue to update VoAGI; if you have any questions or suggestions, please contact us!