合同会社小村ソフト
第 6 章

2 変数ソルバを読む

頂点候補の列挙と可行判定を、外部ライブラリを使わない最小実装で追う。

2 変数ソルバが境界線の交点を列挙して最適点を選ぶ流れ

2 変数なら、境界線の交点を全部調べて、可行な頂点だけを採点できます。

2 変数では「頂点候補を全部見る」ができる

第 4 章で見たように、最適解は通常は頂点に現れます。そこで 2 変数では、次の手順でそのまま解けます。

  1. 境界線 x = 0y = 02x + y = 16x + 2y = 14 を並べる。
  2. 2 本ずつの交点を計算する。
  3. すべての制約を満たす点だけを残す。
  4. 目的関数を評価して最大の点を選ぶ。

この方法は一般の大規模 LP には向きませんが、入門には最適です。

手計算とコードの対応表

これから読む Python コードは、前章までの手計算をそのまま自動化したものです。各ステップの対応関係を先に見ておきます。

手計算でやったこと対応する関数 / コード片
2 本の境界線の交点を連立で求める第 3 章intersection(l1, l2)
点が制約をすべて満たすか代入で確認する第 3 章feasible(point, constraints)
境界線 4 本から 2 本ずつ選ぶ全組合せを並べる第 6 章で新規solve_lp_2d 内の二重ループ
各頂点で目的値 z = c_x x + c_y y を採点する第 4 章max(vertices, key=...)

つまりコードは、第 3 章と第 4 章の作業を「2 本の選び方」のループで包んだだけです。難しいのは Python の文法ではなく、「手計算をどう機械的に並べたか」を読み取る部分です。

境界線をなぜ辞書(dict)で表すのか

このコードでは、境界線を {"kind": "x0"}{"kind": "y0"}{"kind": "line", "a": ..., "b": ..., "c": ...} という辞書で表現します。クラスを作らずに辞書で済ませている理由は次の通りです。

  • 境界線は 3 種類だけ: x = 0(y 軸)、y = 0(x 軸)、一般式 ax + by = c の 3 種類しか登場しません。kind という文字列で見分ければ十分です。
  • 外部依存ゼロを保つため: dataclasses や独自クラスを使わず、組み込み型だけで完結させると、コピペで動かして読み解くのが楽になります。
  • 場合分けで読みやすい: 「軸との交点」と「一般 2 直線の交点」では計算式が違うため、kind で分岐するのが結果的に式に近い表現になります。

大規模なコードベースなら @dataclasstyping.Literal を使うべき場面ですが、教材では場合分けが目に見えるほうが学びやすいのでこの形を採用しています。

最小実装の Python

外部ライブラリを使わず、境界線の交点と可行判定だけで解く最小実装は次のようになります。

EPS = 1e-9  # floating-point tolerance


def intersection(l1, l2):
    kind1 = l1["kind"]
    kind2 = l2["kind"]

    if kind1 == "x0" and kind2 == "y0":
        return (0.0, 0.0)
    if kind1 == "y0" and kind2 == "x0":
        return (0.0, 0.0)

    if kind1 == "x0" and kind2 == "line":
        return (0.0, l2["c"] / l2["b"])
    if kind1 == "line" and kind2 == "x0":
        return (0.0, l1["c"] / l1["b"])

    if kind1 == "y0" and kind2 == "line":
        return (l2["c"] / l2["a"], 0.0)
    if kind1 == "line" and kind2 == "y0":
        return (l1["c"] / l1["a"], 0.0)

    det = l1["a"] * l2["b"] - l2["a"] * l1["b"]
    if abs(det) < EPS:
        return None  # parallel lines: no unique intersection

    x = (l1["c"] * l2["b"] - l2["c"] * l1["b"]) / det
    y = (l1["a"] * l2["c"] - l2["a"] * l1["c"]) / det
    return (x, y)


def feasible(point, constraints):
    x, y = point
    if x < -EPS or y < -EPS:
        return False  # non-negativity x, y >= 0
    return all(a * x + b * y <= c + EPS for a, b, c in constraints)


def solve_lp_2d(cx, cy, constraints):
    if len(constraints) != 2:
        raise ValueError(
            "this 2D solver expects exactly 2 constraints; got "
            + str(len(constraints))
        )

    boundaries = [
        {"kind": "x0"},
        {"kind": "y0"},
        {"kind": "line", "a": constraints[0][0], "b": constraints[0][1], "c": constraints[0][2]},
        {"kind": "line", "a": constraints[1][0], "b": constraints[1][1], "c": constraints[1][2]},
    ]

    candidates = []
    for i in range(len(boundaries)):
        for j in range(i + 1, len(boundaries)):
            p = intersection(boundaries[i], boundaries[j])
            if p is None:
                continue
            if feasible(p, constraints):
                candidates.append((round(p[0], 9), round(p[1], 9)))

    vertices = sorted(set(candidates))
    if not vertices:
        raise ValueError("feasible region is empty")

    best = max(vertices, key=lambda p: cx * p[0] + cy * p[1])
    best_value = cx * best[0] + cy * best[1]
    return best, best_value, vertices

前提が崩れて可行領域が空なら、そのまま ValueError を投げます。2 変数の教材としては、暗黙に処理するより明示的に例外を投げる方が学習時に追跡しやすくなります。

コードの読みどころ

EPS = 1e-9 という閾値の意味

浮動小数点演算では、本来 0 になるべき値が 1e-15 程度ずれることがあります。EPS = 1e-9 は、「これより小さい差は誤差とみなす」という許容範囲です。

  • abs(det) < EPS: 行列式がほぼ 0、つまり 2 直線が平行で交点なしとみなす(1 / 0 によるゼロ除算を避ける)。
  • x < -EPS: 非負制約 x ≥ 0 をわずかに負へ越える誤差は許す(厳密に x < 0 で切ると、本来 0 のはずの頂点が落ちることがある)。
  • a * x + b * y <= c + EPS: 制約も同じく、計算誤差で境界線上の点を弾かないように緩める。

値が 1e-9 なのは「典型的な計算で生じる丸め誤差より十分大きく、実問題のスケールよりは十分小さい」一般的な目安です。係数が 10⁶ 以上のスケールになるなら 1e-6 程度に上げる、といった具合に、扱う数値の桁に応じて調整します。

det = 0(平行線)はなぜスキップするか

行列式 det = a₁b₂ − a₂b₁ が 0 のとき、2 つの境界線は平行(または完全一致)です。平行な 2 直線には交点がない(または無限に多い)ので、頂点候補に加えようがありません。None を返してスキップし、ループ側で if p is None: continue で読み飛ばします。

intersection の場合分けが多い理由

境界線は 3 種類(x0y0line)しかないので、組合せは原理的に 3 × 3 = 9 通り、対称性を考えても 6 通りで、それを愚直に列挙しています。2 変数では境界線が 4 本しかないため、この愚直な列挙でも全体は 4C2 = 6 回しか呼ばれません。一般の n 変数や複数制約の場合はこの方針では破綻するので、教材として「2 変数だから許される単純化」だと意識して読んでください。

sorted(set(candidates)) で重複除去する理由

round(p[0], 9) で丸めてから set にかけているのは、同じ頂点が浮動小数点誤差でわずかに違う値として複数回現れることがあるからです。たとえば (0.0, 0.0)x0 ∩ y0x0 ∩ line(仕上工程の y 切片が偶然 0 のとき)の両方で出てくることがあります。9 桁で丸めれば実用上の差は埋まり、set で重複が消えます。

制約数のバリデーション

len(constraints) != 2 をチェックしているのは、この実装が境界線 4 本を前提にしているからです。3 本以上渡されると boundaries に取りこぼしが出て黙って間違った答えを返してしまうため、入口で例外を投げて気づけるようにしています。

非負制約だけ別扱いする理由

feasible 内では非負制約 x, y ≥ 0constraints 配列に入れず、専用のチェックで処理しています。これは、非負制約は他の制約と扱いが違うわけではないのですが、コード上は x = 0y = 0 を境界線リストにすでに入れているため、constraints 配列には工程制約だけを渡す設計にしている、という整理です。読みやすさのために重複を避けただけで、理論的には他の制約と同じ位置づけです。

理解チェック — 実装の各ステップを読む

実装の各行が、これまでの手計算とどうつながっているかを数で確認します。

Q1. 制約 2x + y = 16x + 2y = 14 を連立するとき、行列式 det = a₁b₂ − a₂b₁ はいくつでしょうか。

Q2. 境界線 x=0y=02x+y=16x+2y=14 の 4 本から 2 本ずつ選ぶ組み合わせは全部でいくつですか。

Q3. 6 個の交点候補のうち、可行判定を通って最終的に残る頂点はいくつですか。

Q4. もし可行判定をせずに点 (0, 16) まで候補に入れてしまうと、z = 3x + 5y の値はいくつになりますか。

万円

第 6 章のまとめ

  • 2 変数では、境界線の交点を全部列挙しても十分小さい。
  • 交点計算のあとに、必ず可行判定で外側の点を落とす。
  • 最終的には、残った頂点で目的関数を比較するだけでよい。