Day 131-Q1 — ネットワーク単体法（輸送問題のMODI法）

2026-08-23 赤色 Master / Phase 8+ ★★★★★★★★★ 北西隅法+MODI法(stepping stone)

問題

$N$個の供給地と$M$個の需要地がある。供給地$i$の供給量は$s_i$、需要地$j$の需要量は$d_j$であり、供給地$i$から需要地$j$へ荷物を1単位輸送するコストは$c_{ij}$である。$\sum_i s_i = \sum_j d_j$が保証されている（均衡輸送問題）。すべての需要をちょうど満たすように輸送量$x_{ij}\ge 0$を決めるとき、総輸送コスト$\sum_{i,j} c_{ij}x_{ij}$の最小値を求めよ。

入力形式

N M
s_1 s_2 ... s_N
d_1 d_2 ... d_M
c_11 c_12 ... c_1M
c_21 c_22 ... c_2M
...
c_N1 c_N2 ... c_NM

制約

$1 \le N, M \le 50$
$1 \le s_i, d_j \le 10^6$
$\sum_i s_i = \sum_j d_j$
$0 \le c_{ij} \le 10^6$

入出力例

入力例1

3 4
7 9 18
5 8 7 14
19 30 50 10
70 30 40 60
40 8 70 20

出力例1

743

概念図

供給地 → 需要地 の二部ネットワーク（基底セルが木を形成） S1 S2 S3 D1 D2 D3 D4 赤: 基底セル(N+M-1本, 木構造) / 灰点線: 非基底セル(被約費用を判定)

ヒント（段階的開示）

ヒント1（方向性）

これは輸送問題という特殊構造を持つ最小費用流であり、一般のMCMFでも解けるが、二部グラフに特化した輸送単体法（Transportation Simplex＝ネットワーク単体法の特化版）を使うと、単体法の各ステップが「表の上のセルの出入り」として直感的に実行できる。

ヒント2（アプローチ）

①北西隅法で初期基底可能解（N+M-1個の基底セル）を作る。②基底セルだけから双対変数u_i,v_jを求める（基底で c_ij=u_i+v_j）。③非基底セルの被約費用 c_ij-(u_i+v_j) が負なら、そのセルを基底に入れ、既存基底との唯一の閉路に沿って輸送量を交互に増減させる（stepping stone）。すべて非負になったら最適。

ヒント3（誘導） 
def find_loop(cells_set, start):
    cells = set(cells_set); cells.add(start)
    changed = True
    while changed:
        changed = False
        for c in list(cells):
            row_cnt = sum(1 for cc in cells if cc[0] == c[0])
            col_cnt = sum(1 for cc in cells if cc[1] == c[1])
            if row_cnt < 2 or col_cnt < 2:
                cells.discard(c); changed = True
    # cells に残ったものが唯一の閉路

模範解答 (Python)

import sys

def solve_transportation(supply, demand, cost):
    m, n = len(supply), len(demand)
    x = [[0] * n for _ in range(m)]
    basis = set()

    s, d = supply[:], demand[:]
    i = j = 0
    while i < m and j < n:
        amt = min(s[i], d[j])
        x[i][j] = amt
        basis.add((i, j))
        s[i] -= amt
        d[j] -= amt
        if s[i] == 0 and d[j] == 0:
            if i < m - 1:
                i += 1
            elif j < n - 1:
                j += 1
            else:
                break
        elif s[i] == 0:
            i += 1
        else:
            j += 1

    def get_uv():
        u = [None] * m
        v = [None] * n
        u[0] = 0
        changed = True
        while changed:
            changed = False
            for (bi, bj) in basis:
                if u[bi] is not None and v[bj] is None:
                    v[bj] = cost[bi][bj] - u[bi]
                    changed = True
                elif v[bj] is not None and u[bi] is None:
                    u[bi] = cost[bi][bj] - v[bj]
                    changed = True
        return u, v

    def find_loop(start):
        cells = set(basis)
        cells.add(start)
        changed = True
        while changed:
            changed = False
            for c in list(cells):
                row_cnt = sum(1 for cc in cells if cc[0] == c[0])
                col_cnt = sum(1 for cc in cells if cc[1] == c[1])
                if row_cnt < 2 or col_cnt < 2:
                    cells.discard(c)
                    changed = True
        order = [start]
        used = {start}
        cur = start
        horizontal = True
        while True:
            if horizontal:
                nxts = [c for c in cells if c[0] == cur[0] and c not in used]
            else:
                nxts = [c for c in cells if c[1] == cur[1] and c not in used]
            if not nxts:
                break
            cur = nxts[0]
            order.append(cur)
            used.add(cur)
            horizontal = not horizontal
        return order

    for _ in range(2000):
        u, v = get_uv()
        enter, best = None, 0
        for a in range(m):
            for b in range(n):
                if (a, b) not in basis:
                    rc = cost[a][b] - (u[a] + v[b])
                    if rc < best:
                        best, enter = rc, (a, b)
        if enter is None:
            break  # 最適
        loop = find_loop(enter)
        minus_cells = loop[1::2]
        theta = min(x[a][b] for (a, b) in minus_cells)
        leaving = None
        for idx, (a, b) in enumerate(loop):
            if idx % 2 == 0:
                x[a][b] += theta
            else:
                x[a][b] -= theta
                if x[a][b] == 0 and leaving is None:
                    leaving = (a, b)
        basis.add(enter)
        basis.discard(leaving)

    return sum(cost[a][b] * x[a][b] for a in range(m) for b in range(n))


def solve():
    data = sys.stdin.read().split()
    idx = 0
    n = int(data[idx]); idx += 1
    m = int(data[idx]); idx += 1
    supply = [int(data[idx + k]) for k in range(n)]; idx += n
    demand = [int(data[idx + k]) for k in range(m)]; idx += m
    cost = []
    for _ in range(n):
        cost.append([int(data[idx + k]) for k in range(m)])
        idx += m
    print(solve_transportation(supply, demand, cost))


solve()

Step-by-Step 解説

1北西隅法で初期基底可能解
表の左上から順に供給・需要の小さい方を割り当て、消費し尽くしたら次の行/列へ。N+M-1個の基底セル（木構造）が得られる。
2双対変数u,vを求める
基底セルではc_ij=u_i+v_j。基底が木なのでu_1=0と固定すれば残りは一意に伝播で決まる。
3MODI法で最適性判定
非基底セルの被約費用c_ij-(u_i+v_j)が負なら改善余地あり。すべて0以上なら最適（線形計画の相補性定理）。
4閉路に沿ったピボット
被約費用が最小のセルを選び、既存基底との唯一の閉路上で+θ,-θを交互に加える。θ=(-θ側セルの現在値の最小)。最初に0になったセルを基底から外す。
5計算量
基底サイズはN+M-1で固定。ピボット回数は最悪O(NM)程度、各ステップO(NM)で全体O(N²M²)程度。N,M≤50なら十分高速。

よくあるミス

ミス原因正しい書き方
初期基底セル数がN+M-1にならない北西隅法の終端で供給・需要が同時に0になる分岐漏れ同時に尽きたら最後の行/列でなければ片方だけ進める特別処理を入れる
双対変数の基準点を毎回変えるu[0]=0を固定せず不定形の連立方程式を解こうとする根を1つ固定すれば基底の木構造上でBFS的に一意に伝播する
閉路上の符号を逆にするentering cellから交互に+,-を割り当てる順序ミス0番目(entering)を+θ、奇数番目を-θと固定してからθを決める
不均衡輸送問題を想定してしまうダミー供給地/需要地の追加を忘れる本問は均衡が保証されるが、実務では差分を吸収するコスト0のダミー行/列を追加する

次のステップ

  • 発展: 容量制約付き輸送問題や一般グラフ上の最小費用流に拡張したネットワーク単体法（スパニングツリー基底+ポテンシャル更新+木を跨ぐ辺のピボット）
  • 次回予告: 次回のテーマは今後のセッションで選定

自己評価

自分の回答

気づき・メモ