prob.writeLP("my_model.lp")

import numpy as np
from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt

# 問題空間の描画
def plot_vrp_nodes(coords):
    plt.figure(figsize=(8, 8))
    
    # 顧客ノード（1以降）を青で描画
    plt.scatter(coords[1:, 0], coords[1:, 1], c="tab:blue", label="Customers")

    # デポ（0番）を赤い星で描画
    plt.scatter(coords[0, 0], coords[0, 1], c="red", s=100, marker="*", label="Depot") # type: ignore

    # ラベル（ノード番号）を各点に表示
    for i, (x, y) in enumerate(coords):
        plt.text(x + 1, y + 1, str(i), fontsize=8)
    
    plt.title(f"VRP Nodes (1 Depot + {len(coords) - 1} Customers)")
    plt.xlabel("X Coordinate")
    plt.ylabel("Y Coordinate")
    plt.legend()
    plt.grid(True)
    plt.axis("equal")
    plt.tight_layout()
    plt.show()


np.random.seed(0)  # 再現性のためのシード固定

# --- パラメータ設定 ---
n_customers = 10         # 顧客の数
n_total = n_customers + 1  # 総ノード数（デポ1 + 顧客）
coord_range = (0, 100)     # 座標の範囲

# --- ノード座標の生成 ---
# 座標: 配列 [ [x0, y0], [x1, y1], ... ]  (0番がデポ)
coords = np.random.uniform(low=coord_range[0], high=coord_range[1], size=(n_total, 2))

# --- ユークリッド距離行列の計算 ---
dist_matrix = cdist(coords, coords)  # shape = (n_total, n_total)

# --- チェック ---
print("Depot座標:", coords[0])

# --- 問題空間の描画 ---
plot_vrp_nodes(coords)

Depot座標: [54.88135039 71.51893664]

import numpy as np
import pulp
import itertools
from tqdm import tqdm

# ノード集合
nodes = list(range(n_total))  # 0: デポ, 1〜n_customers: 顧客
n = n_customers # 顧客数
m = 2  # 車両台数

# 距離行列（対称）
distance_matrix = dist_matrix


# 最適化問題の定義
prob = pulp.LpProblem("VRP_with_SEC", pulp.LpMinimize)

# バイナリ変数 x_{ij}
x = {(i, j): pulp.LpVariable(f"x_{i}_{j}", cat=pulp.LpBinary)
     for i in nodes for j in nodes if i != j}

# 目的関数：移動距離の合計最小化
prob += pulp.lpSum(distance_matrix[i][j] * x[i, j] for (i, j) in x)

# 出発と帰着の台数制約
prob += pulp.lpSum(x[0, j] for j in nodes if j != 0) == m
prob += pulp.lpSum(x[i, 0] for i in nodes if i != 0) == m

# 流入と流出の制約（各顧客ノード）
for k in nodes:
    if k == 0:
        continue
    prob += pulp.lpSum(x[i, k] for i in nodes if i != k) == 1
    prob += pulp.lpSum(x[k, j] for j in nodes if j != k) == 1


# 顧客ノード集合（0は除く）
customers = [i for i in nodes if i != 0]

# 全ての部分集合 S に対して SEC を追加する(顧客数が小さいので全列挙できる)
for r in tqdm(range(2, len(customers) + 1)):
    for S in itertools.combinations(customers, r):
        prob += pulp.lpSum(x[i, j] for i in S for j in S if i != j) <= len(S) - 1

# 解の実行と表示
from pulp import PULP_CBC_CMD

# 求解
solver = PULP_CBC_CMD(msg=True)
prob.solve(solver)

# 結果表示
for (i, j) in x:
    if x[i, j].varValue == 1:
        print(f"Vehicle travels from {i} to {j}")

100%|██████████| 9/9 [00:00<00:00, 670.37it/s]

Welcome to the CBC MILP Solver 
Version: 2.10.10 
Build Date: Sep 26 2023 

command line - /home/kakik/.pyenv/versions/3.14.0/lib/python3.14/site-packages/pulp/apis/../solverdir/cbc/linux/arm64/cbc /tmp/eefda76438f74f1f94276d363e0833cd-pulp.mps -timeMode elapsed -branch -printingOptions all -solution /tmp/eefda76438f74f1f94276d363e0833cd-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 1040 COLUMNS
At line 24631 RHS
At line 25667 BOUNDS
At line 25778 ENDATA
Problem MODEL has 1035 rows, 110 columns and 23260 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 381.741 - 0.00 seconds
Cgl0004I processed model has 1035 rows, 110 columns (110 integer (110 of which binary)) and 23260 elements
Cbc0038I Initial state - 0 integers unsatisfied sum - 0
Cbc0038I Solution found of 381.741
Cbc0038I Before mini branch and bound, 110 integers at bound fixed and 0 continuous
Cbc0038I Mini branch and bound did not improve solution (0.02 seconds)
Cbc0038I After 0.02 seconds - Feasibility pump exiting with objective of 381.741 - took 0.00 seconds
Cbc0012I Integer solution of 381.74063 found by feasibility pump after 0 iterations and 0 nodes (0.02 seconds)
Cbc0001I Search completed - best objective 381.74063484596, took 0 iterations and 0 nodes (0.02 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from 381.741 to 381.741
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                381.74063485
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.01
Time (Wallclock seconds):       0.03

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.02   (Wallclock seconds):       0.03

Vehicle travels from 0 to 1
Vehicle travels from 0 to 6
Vehicle travels from 1 to 5
Vehicle travels from 2 to 0
Vehicle travels from 3 to 8
Vehicle travels from 4 to 10
Vehicle travels from 5 to 4
Vehicle travels from 6 to 3
Vehicle travels from 7 to 2
Vehicle travels from 8 to 7
Vehicle travels from 9 to 0
Vehicle travels from 10 to 9

def plot_vrp_solution(coords, x):
    """
    coords : ndarray of shape (n_nodes, 2)
        ノードの座標。0番がデポ。
    x : dict {(i, j): pulp.LpVariable}
        解を持つバイナリ変数。x[i, j].varValue == 1 ならルート (i→j) を通る。
    """
    plt.figure(figsize=(8, 8))

    # ノードを描画（デポ：赤星、顧客：青点）
    plt.scatter(coords[1:, 0], coords[1:, 1], c="tab:blue", label="Customers")
    plt.scatter(coords[0, 0], coords[0, 1], c="red", s=100, marker="*", label="Depot (0)")

    # ラベルを表示
    for i, (x_, y_) in enumerate(coords):
        plt.text(x_ + 1, y_ + 1, str(i), fontsize=8)

    # 解からルートを描画
    for (i, j), var in x.items():
        if hasattr(var, 'varValue') and var.varValue > 0.5:
            xi, yi = coords[i]
            xj, yj = coords[j]
            plt.plot([xi, xj], [yi, yj], 'k-', alpha=0.5)

    plt.title("VRP Solution")
    plt.xlabel("X Coordinate")
    plt.ylabel("Y Coordinate")
    plt.legend()
    plt.grid(True)
    plt.axis("equal")
    plt.tight_layout()
    plt.show()

# 最適解を描画（xは {(i, j): LpVariable} の辞書）
plot_vrp_solution(coords, x)

import networkx as nx

# ダミー道路網の読み込み(road_network.ipynb と同じ手順)
edgs = np.loadtxt('../road-network/links_ud.csv', delimiter=',')
G = nx.Graph()
G.add_weighted_edges_from([(int(e[0]), int(e[1]), e[2]) for e in edgs])
pos_nd = {}
for e in edgs:
    pos_nd[int(e[0])] = [e[3], e[4]]
    pos_nd[int(e[1])] = [e[5], e[6]]

# 地点番号と道路網ノードIDの対応表(0番がデポ)
spots = [15, 1, 6, 13, 24, 25, 30]
n2 = len(spots)
m2 = 2  # 車両台数

# 最短経路距離による距離行列
dist2 = np.zeros((n2, n2))
for a in range(n2):
    for b in range(n2):
        if a != b:
            dist2[a][b] = nx.shortest_path_length(G, spots[a], spots[b], weight='weight')

print('地点間の最短経路距離 [m](行と列の並びは spots の順):')
print(np.round(dist2, 1))

地点間の最短経路距離 [m](行と列の並びは spots の順):
[[   0.  1336.8 1496.2  816.1 1180.  1315.7 1542.9]
 [1336.8    0.  1787.6 1484.3 2306.  2034.4 2668.9]
 [1496.2 1787.6    0.  2312.3 1046.9 2811.9 1409.8]
 [ 816.1 1484.3 2312.3    0.  1996.1  718.  2359. ]
 [1180.  2306.  1046.9 1996.1    0.  1991.4  362.9]
 [1315.7 2034.4 2811.9  718.  1991.4    0.  1882.7]
 [1542.9 2668.9 1409.8 2359.   362.9 1882.7    0. ]]

nodes2 = list(range(n2))
prob2 = pulp.LpProblem('VRP_road_network', pulp.LpMinimize)
x2 = {(i, j): pulp.LpVariable(f'x2_{i}_{j}', cat=pulp.LpBinary)
      for i in nodes2 for j in nodes2 if i != j}

prob2 += pulp.lpSum(dist2[i][j] * x2[i, j] for (i, j) in x2)

prob2 += pulp.lpSum(x2[0, j] for j in nodes2 if j != 0) == m2
prob2 += pulp.lpSum(x2[i, 0] for i in nodes2 if i != 0) == m2
for k in nodes2:
    if k == 0:
        continue
    prob2 += pulp.lpSum(x2[i, k] for i in nodes2 if i != k) == 1
    prob2 += pulp.lpSum(x2[k, j] for j in nodes2 if j != k) == 1
customers2 = [i for i in nodes2 if i != 0]
for r in range(2, len(customers2) + 1):
    for S in itertools.combinations(customers2, r):
        prob2 += pulp.lpSum(x2[i, j] for i in S for j in S if i != j) <= len(S) - 1

prob2.solve(PULP_CBC_CMD(msg=0))
print('求解結果:', pulp.LpStatus[prob2.status])
print('総移動距離:', round(pulp.value(prob2.objective), 1), '[m]')

# ルートの復元(デポから出るアークをたどり、ノードIDに変換して表示)
succ = {i: j for (i, j) in x2 if i != 0 and x2[i, j].varValue > 0.5}
routes = []
for s in [j for j in customers2 if x2[0, j].varValue > 0.5]:
    route = [0, s]
    while route[-1] != 0:
        route.append(succ[route[-1]])
    routes.append(route)
for rt in routes:
    print('ルート:', ' -> '.join(str(spots[k]) for k in rt))

求解結果: Optimal
総移動距離: 8926.9 [m]
ルート: 15 -> 1 -> 6 -> 24 -> 30 -> 15
ルート: 15 -> 25 -> 13 -> 15

best_val = None
for mask in range(1, 2 ** len(customers2) - 1):
    g1 = [customers2[k] for k in range(len(customers2)) if mask >> k & 1]
    g2 = [c for c in customers2 if c not in g1]
    total = 0
    for g in (g1, g2):
        total += min(dist2[0][perm[0]]
                     + sum(dist2[perm[i]][perm[i + 1]] for i in range(len(perm) - 1))
                     + dist2[perm[-1]][0]
                     for perm in itertools.permutations(g))
    if best_val is None or total < best_val:
        best_val = total

print('総当たりによる最小総移動距離:', round(best_val, 1), '[m]')
print('数理モデルの解と一致:', round(best_val, 1) == round(pulp.value(prob2.objective), 1))

総当たりによる最小総移動距離: 8926.9 [m]
数理モデルの解と一致: True

# 可視化: 各車両のルートを、道路上の実際の最短経路に沿って描く
route_colors = ['red', 'tab:blue']
plt.figure(figsize=(8, 6))
nx.draw_networkx_edges(G, pos_nd, edge_color='lightgray')
for color, route in zip(route_colors, routes):
    for a, b in zip(route[:-1], route[1:]):
        sp = nx.shortest_path(G, spots[a], spots[b], weight='weight')
        nx.draw_networkx_edges(G, pos_nd, edgelist=list(zip(sp[:-1], sp[1:])),
                               edge_color=color, width=2.5)
nx.draw_networkx_nodes(G, pos_nd, node_size=170, node_color='white', edgecolors='gray')
nx.draw_networkx_nodes(G, pos_nd, nodelist=[spots[i] for i in customers2], node_size=260,
                       node_color='lightblue', edgecolors='black')
nx.draw_networkx_nodes(G, pos_nd, nodelist=[spots[0]], node_shape='*', node_size=600,
                       node_color='red', edgecolors='black')
nx.draw_networkx_labels(G, pos_nd, font_size=7)
plt.axis('equal')
plt.axis('off')
plt.tight_layout()
plt.show()

PuLPのコード例	LPファイル上の出力例
`prob += lpSum(c[i,j] * x[i,j] for ...)`	`obj: ...` の式
`prob += lpSum(x[i,j] for ...) == 1`	制約条件 `c1: ... = 1`
`LpVariable(..., cat=LpBinary)`	`Binary` 宣言

Vehicle Routing Problem (VRP) 入門¶

1. はじめに¶

2. 前提知識と学習目標¶

3. 問題設定と定式化¶

3.1 問題の具体例¶

3.2 数学的な定式化¶

3.2.1 定数および決定変数の定義¶

3.2.2 数理モデル¶

3.2.3 目的関数の意味¶

3.2.4 制約条件の意味¶

4. LPファイルとの対応関係¶

4.1 LPファイルの基本構造¶

4.2 PuLPのコードとの対応¶

5. Python + PuLP による実装¶

6. 道路ネットワーク上のVRP¶

7. 注意事項¶

7.1 部分巡回路除去制約（サブツアー制約）について¶

7.2 定式化が理解できないときは¶

8. 卒業研究概要への手引き¶

1章はじめに(文章量の比率: 15%)¶

2章 VRP(文章量の比率: 30%)¶

3章実験(文章量の比率: 40%)¶

3.1節概要¶

3.2節結果と考察¶

4章おわりに(文章量の比率: 15%)¶

参考文献¶

Vehicle Routing Problem (VRP) 入門¶

1. はじめに¶

2. 前提知識と学習目標¶

3. 問題設定と定式化¶

3.1 問題の具体例¶

3.2 数学的な定式化¶

3.2.1 定数および決定変数の定義¶

3.2.2 数理モデル¶

3.2.3 目的関数の意味¶

3.2.4 制約条件の意味¶

4. LPファイルとの対応関係¶

4.1 LPファイルの基本構造¶

4.2 PuLPのコードとの対応¶

5. Python + PuLP による実装¶

6. 道路ネットワーク上のVRP¶

7. 注意事項¶

7.1 部分巡回路除去制約（サブツアー制約）について¶

7.2 定式化が理解できないときは¶

8. 卒業研究概要への手引き¶

1章 はじめに(文章量の比率: 15%)¶

2章 VRP(文章量の比率: 30%)¶

3章 実験(文章量の比率: 40%)¶

3.1節 概要¶

3.2節 結果と考察¶

4章 おわりに(文章量の比率: 15%)¶

参考文献¶

1章はじめに(文章量の比率: 15%)¶

3章実験(文章量の比率: 40%)¶

3.1節概要¶

3.2節結果と考察¶

4章おわりに(文章量の比率: 15%)¶