import itertools
import pulp

d = [0, 5, 10, 15, 20]  # ノード1〜5の位置 [km](d[i-1] がノードi の位置)
K = [2, 3, 4]           # 施設の候補地
N = [[1, 2, 3, 4, 5], [1, 2, 3], [2, 3, 4, 5], [3, 4, 5]]  # 経路の通過ノード列
f = [100, 60, 50, 40]                                      # 経路の交通量 [台/日]
J = list(range(1, len(N) + 1))
R = 12  # 航続距離 [km]
p = 2   # 施設数


def refuelable(Nj, combo, R):
    """経路Nj の往復が、施設の組合せcombo で補給可能かを判定する"""
    org = d[Nj[0] - 1]
    L = d[Nj[-1] - 1] - org
    pts = [0] + sorted(d[k - 1] - org for k in combo)
    for a, b in zip(pts[:-1], pts[1:]):
        if b - a > R:                 # 条件1、2
            return False
    return 2 * (L - pts[-1]) <= R     # 条件3


# 各経路について補給可能な組合せ H_j を列挙する
H = []
for j in J:
    on_path = [k for k in K if k in N[j - 1]]
    combos = [h for r in range(1, len(on_path) + 1)
              for h in itertools.combinations(on_path, r)
              if refuelable(N[j - 1], h, R)]
    H.append(combos)
    print(f'経路{j}: 補給可能な組合せ {combos}')

経路1: 補給可能な組合せ [(2, 4), (3, 4), (2, 3, 4)]
経路2: 補給可能な組合せ [(2,), (3,), (2, 3)]
経路3: 補給可能な組合せ [(4,), (2, 4), (3, 4), (2, 3, 4)]
経路4: 補給可能な組合せ [(4,), (3, 4)]

all_combos = sorted({h for combos in H for h in combos})

prob = pulp.LpProblem('FRLM', pulp.LpMaximize)
x = pulp.LpVariable.dicts('x', K, cat='Binary')
v = pulp.LpVariable.dicts('v', range(len(all_combos)), cat='Binary')
y = pulp.LpVariable.dicts('y', J, cat='Binary')

prob += pulp.lpSum(f[j - 1] * y[j] for j in J)
for j in J:
    prob += y[j] <= pulp.lpSum(v[all_combos.index(h)] for h in H[j - 1])
for i, h in enumerate(all_combos):
    for k in h:
        prob += v[i] <= x[k]
prob += pulp.lpSum(x[k] for k in K) == p

prob.solve(pulp.PULP_CBC_CMD(msg=0))
print('求解結果:', pulp.LpStatus[prob.status])
print('捕捉交通量:', int(pulp.value(prob.objective)), '/', sum(f))
print('開設ノード:', [k for k in K if x[k].value() == 1])
print('補給可能な経路:', [j for j in J if y[j].value() == 1])

求解結果: Optimal
捕捉交通量: 250 / 250
開設ノード: [2, 4]
補給可能な経路: [1, 2, 3, 4]

def solve_fclm(p):
    prob = pulp.LpProblem('FCLM', pulp.LpMaximize)
    x = pulp.LpVariable.dicts('x', K, cat='Binary')
    y = pulp.LpVariable.dicts('y', J, cat='Binary')
    prob += pulp.lpSum(f[j - 1] * y[j] for j in J)
    for j in J:
        prob += y[j] <= pulp.lpSum(x[k] for k in K if k in N[j - 1])
    prob += pulp.lpSum(x[k] for k in K) == p
    prob.solve(pulp.PULP_CBC_CMD(msg=0))
    return int(pulp.value(prob.objective)), [k for k in K if x[k].value() == 1]


def solve_frlm(p, R):
    H = []
    for j in J:
        on_path = [k for k in K if k in N[j - 1]]
        H.append([h for r in range(1, len(on_path) + 1)
                  for h in itertools.combinations(on_path, r)
                  if refuelable(N[j - 1], h, R)])
    all_combos = sorted({h for combos in H for h in combos})
    prob = pulp.LpProblem('FRLM', pulp.LpMaximize)
    x = pulp.LpVariable.dicts('x', K, cat='Binary')
    v = pulp.LpVariable.dicts('v', range(len(all_combos)), cat='Binary')
    y = pulp.LpVariable.dicts('y', J, cat='Binary')
    prob += pulp.lpSum(f[j - 1] * y[j] for j in J)
    for j in J:
        prob += y[j] <= pulp.lpSum(v[all_combos.index(h)] for h in H[j - 1])
    for i, h in enumerate(all_combos):
        for k in h:
            prob += v[i] <= x[k]
    prob += pulp.lpSum(x[k] for k in K) == p
    prob.solve(pulp.PULP_CBC_CMD(msg=0))
    return int(pulp.value(prob.objective)), [k for k in K if x[k].value() == 1]


print('p | FCLM(航続距離を無視) | FRLM(R=12)')
for p_test in [1, 2, 3]:
    z1, s1 = solve_fclm(p_test)
    z2, s2 = solve_frlm(p_test, R)
    print(f'{p_test} | {z1} {s1} | {z2} {s2}')

p | FCLM(航続距離を無視) | FRLM(R=12)
1 | 250 [3] | 90 [4]
2 | 250 [2, 3] | 250 [2, 4]
3 | 250 [2, 3, 4] | 250 [2, 3, 4]

print('R | 捕捉交通量(p=2) | 開設ノード')
for R_test in [8, 10, 12, 16, 20]:
    z, s = solve_frlm(2, R_test)
    print(f'{R_test} | {z} | {s}')

R | 捕捉交通量(p=2) | 開設ノード
8 | 60 | [2, 3]
10 | 250 | [2, 4]
12 | 250 | [2, 4]
16 | 250 | [2, 4]
20 | 250 | [2, 3]

import numpy as np
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 = {}
for e in edgs:
    pos[int(e[0])] = [e[3], e[4]]
    pos[int(e[1])] = [e[5], e[6]]

# 経路データの読み込み
N = []
with open('../road-network/N.csv', 'r') as fp:
    for row in fp.readlines():
        N.append([int(v) for v in row.strip().split(',') if len(v) > 0])
f = [int(v) for v in np.loadtxt('../road-network/f.csv')]
J = list(range(1, len(N) + 1))
I = sorted(set(sum(N, [])))
R_nw = 1600.0  # 航続距離 [m]


def cum_dist(Nj):
    """経路のノード列に沿った累積距離"""
    dd = [0.0]
    for u, v in zip(Nj[:-1], Nj[1:]):
        dd.append(dd[-1] + G[u][v]['weight'])
    return dd


def refuelable_path(Nj, combo, R):
    """経路Nj の往復が、施設の組合せcombo で補給可能かを判定する"""
    dd = cum_dist(Nj)
    L = dd[-1]
    pts = [0] + [dd[i] for i, nd in enumerate(Nj) if nd in combo]
    for a, b in zip(pts[:-1], pts[1:]):
        if b - a > R:                 # 条件1、2
            return False
    return 2 * (L - pts[-1]) <= R     # 条件3


# 各経路について補給可能な組合せ H_j を列挙する(候補は経路上の全ノード)
H_nw = []
for j in J:
    on_path = N[j - 1]
    combos = [h for r in range(len(on_path) + 1)
              for h in itertools.combinations(on_path, r)
              if refuelable_path(N[j - 1], h, R_nw)]
    H_nw.append(combos)
    L = cum_dist(N[j - 1])[-1]
    print(f'経路{j}: 経路長 {L:6.1f} [m], 補給可能な組合せの数 {len(combos)}')

経路1: 経路長 1787.6 [m], 補給可能な組合せの数 54
経路2: 経路長 1882.7 [m], 補給可能な組合せの数 46
経路3: 経路長 2034.4 [m], 補給可能な組合せの数 104
経路4: 経路長 1409.8 [m], 補給可能な組合せの数 28
経路5: 経路長 1292.0 [m], 補給可能な組合せの数 28
経路6: 経路長 1911.7 [m], 補給可能な組合せの数 54
経路7: 経路長 1161.3 [m], 補給可能な組合せの数 14

def solve_frlm_network(p, R):
    H = []
    for j in J:
        on_path = N[j - 1]
        H.append([h for r in range(len(on_path) + 1)
                  for h in itertools.combinations(on_path, r)
                  if refuelable_path(N[j - 1], h, R)])
    all_combos = sorted({h for combos in H for h in combos})
    prob = pulp.LpProblem('FRLM_network', pulp.LpMaximize)
    x = pulp.LpVariable.dicts('x', I, cat='Binary')
    v = pulp.LpVariable.dicts('v', range(len(all_combos)), cat='Binary')
    y = pulp.LpVariable.dicts('y', J, cat='Binary')
    prob += pulp.lpSum(f[j - 1] * y[j] for j in J)
    for j in J:
        prob += y[j] <= pulp.lpSum(v[all_combos.index(h)] for h in H[j - 1])
    for i, h in enumerate(all_combos):
        for k in h:
            prob += v[i] <= x[k]
    prob += pulp.lpSum(x[i] for i in I) == p
    prob.solve(pulp.PULP_CBC_CMD(msg=0))
    z = int(pulp.value(prob.objective))
    chosen = [i for i in I if x[i].value() == 1]
    captured = [j for j in J if y[j].value() == 1]
    return z, chosen, captured


z2, chosen2, captured2 = solve_frlm_network(2, R_nw)
print('p = 2 の捕捉交通量:', z2, '/', sum(f))
print('開設ノード:', chosen2)
print('補給可能な経路:', captured2)

p = 2 の捕捉交通量: 360 / 670
開設ノード: [15, 18]
補給可能な経路: [4, 5, 6, 7]

from itertools import combinations

def frlm_value(open_nodes, R):
    return sum(fj for Nj, fj in zip(N, f)
               if refuelable_path(Nj, set(Nj) & set(open_nodes), R))

best = max(frlm_value(cb, R_nw) for cb in combinations(I, 2))
print('全列挙による p = 2 の最大捕捉交通量:', best)
print('数理モデルの解と一致:', best == z2)

全列挙による p = 2 の最大捕捉交通量: 360
数理モデルの解と一致: True

import matplotlib.pyplot as plt

# 可視化: 赤い実線 = 補給可能な経路、青い破線 = 補給できない経路
plt.figure(figsize=(8, 6))
nx.draw_networkx_edges(G, pos, edge_color='lightgray')
for j in J:
    if j not in captured2:
        nx.draw_networkx_edges(G, pos, edgelist=list(zip(N[j - 1][:-1], N[j - 1][1:])),
                               edge_color='steelblue', style='dashed', width=2)
for j in J:
    if j in captured2:
        nx.draw_networkx_edges(G, pos, edgelist=list(zip(N[j - 1][:-1], N[j - 1][1:])),
                               edge_color='red', width=2)
nx.draw_networkx_nodes(G, pos, node_size=170, node_color='white', edgecolors='gray')
nx.draw_networkx_nodes(G, pos, nodelist=chosen2, node_size=300, node_color='red',
                       edgecolors='black')
nx.draw_networkx_labels(G, pos, font_size=7)
plt.axis('equal')
plt.axis('off')
plt.tight_layout()
plt.show()

for p_test in [1, 2, 3, 4]:
    z, chosen, _ = solve_frlm_network(p_test, R_nw)
    print(f'p = {p_test}: 捕捉交通量 = {z} [台/日], 開設ノード = {chosen}')

print()
print('R [m] | 捕捉交通量(p = 2) | 開設ノード')
for R_test in [1200, 1600, 2000, 2400, 2800]:
    z, chosen, _ = solve_frlm_network(2, R_test)
    print(f'{R_test} | {z} | {chosen}')

p = 1: 捕捉交通量 = 150 [台/日], 開設ノード = [18]
p = 2: 捕捉交通量 = 360 [台/日], 開設ノード = [15, 18]

p = 3: 捕捉交通量 = 480 [台/日], 開設ノード = [4, 15, 18]

p = 4: 捕捉交通量 = 580 [台/日], 開設ノード = [4, 15, 18, 29]

R [m] | 捕捉交通量(p = 2) | 開設ノード
1200 | 300 | [21, 24]
1600 | 360 | [15, 18]

2000 | 360 | [15, 18]
2400 | 460 | [15, 30]

2800 | 520 | [6, 30]

燃料補給施設配置問題(FRLM)入門¶

1. はじめに¶

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

学習目標¶

3. 問題設定と定式化¶

3.1 問題の具体例¶

3.2 補給可能の条件¶

3.3 定式化¶

4. Python + PuLP による実装¶

4.1 データと補給可能判定¶

4.2 モデル構築と求解¶

4.3 FCLMとの比較¶

4.4 感度分析: 航続距離を変える¶

5. 道路ネットワークへの適用¶

6. 発展: 研究紹介¶

7. 参考文献¶

8. 課題¶

9. 注意事項¶

経路 $j$	通過ノード $N_j$	交通量 $f_j$ [台/日]
1	1, 2, 3, 4, 5	100
2	1, 2, 3	60
3	2, 3, 4, 5	50
4	3, 4, 5	40