線形計画による生産計画｜オペレーションズマネジメント

🎓 レベル：標準　|　重要度：A（必須）

📎 前提：ボトルネックとキャパシティ（拘束制約＝ボトルネック・資源キャパシティ）　|　最適化理論：凸最適化の基礎（凸性・双対・KKT。LP はその線形な特殊例）　|　次：スケジューリングとジョブショップ（順序の最適化）

要点（BLUF）

製品ミックス問題：限られた資源（機械時間・原料・人手）を複数製品にどう振り分ければ利益が最大になるか。決定変数を各製品の生産量 $x_i$ 、目的を利益最大 $\max\sum_i c_i x_i$ 、制約を資源 $\sum_i a_{ij}x_i\le b_j$ と $x_i\ge0$ と書けば、これは**線形計画（LP）**そのものです。
実行可能領域（制約をすべて満たす $x$ の集合）は凸多面体。線形の目的関数を動かすと、最適解は必ずその**頂点（端点）**に来ます。だから LP の答えは「どの制約が同時に効いて交わる角か」を選ぶ問題です。
ある制約がぴったり等号で効いている（余り slack $=0$ ）とき、それが拘束制約＝ボトルネック（ボトルネックとキャパシティ）。余っている資源（slack $>0$ ）は律速していません。
シャドープライス（双対価格）＝その資源を1単位増やしたとき増える利益。拘束資源だけが正のシャドープライスを持ち、非拘束資源のシャドープライスは0（余っているものを足しても儲からない＝補完スラック）。これが「どの資源に投資すべきか」を教えます。
LP は凸最適化の線形な特殊例です。凸性・双対・KKT条件の一般理論は凸最適化の基礎に譲り、本ノートは「LP を生産計画にどう使うか」を scipy.optimize.linprog で計算して見せます。

1. 製品ミックス問題を立式する

工場が製品A・製品Bをつくっています。それぞれ1個あたりの利益（貢献利益）は $c_A,c_B$ 。生産には資源（機械加工の時間、組立の時間、原料）が要り、各製品が1個あたり消費する量と、資源の月間上限が決まっています。利益を最大にする生産量の組 $(x_A,x_B)$ を求めたい——これが製品ミックス問題です。

要素は3つだけです。

決定変数 $x_i\ge0$ ：製品 $i$ の生産量（いくつつくるか）。負の生産はないので非負。
目的関数 $\max\ Z=\sum_i c_i x_i$ ：総利益。各製品の利益 × 生産量の和で、 $x$ の1次式。
制約 $\sum_i a_{ij}x_i\le b_j$ ：資源 $j$ の使用量が上限 $b_j$ を超えない。 $a_{ij}$ は製品 $i$ が資源 $j$ を1個あたり消費する量（技術係数）。これも1次式。

目的も制約もすべて1次（線形）なので「線形計画」。EOQ や M/M/1 のように閉形式の公式は出ませんが、その代わり変数や制約がいくつあってもソルバが頂点を探し当てます。

flowchart LR
  X["決定変数<br/>各製品の生産量 x_i >= 0"] --> OBJ["目的<br/>利益最大 max Σ c_i x_i"]
  X --> CON["資源制約<br/>Σ a_ij x_i <= b_j"]
  OBJ --> LP["線形計画 LP"]
  CON --> LP
  LP --> V["最適は実行可能領域の頂点<br/>拘束制約 slack=0 がボトルネック"]

2. なぜ最適は「頂点」なのか・シャドープライスとは

実行可能領域は、各制約 $\sum_i a_{ij}x_i\le b_j$ が定める半空間の共通部分です。半空間は凸で、凸集合の共通部分も凸なので、実行可能領域は凸多面体（角のある多角形・多面体）になります。

目的関数 $Z=c^\top x$ は線形なので、 $Z=\text{一定}$ の等高線（等利益線）は平行な直線（超平面）の束。これを利益が増える向きへ平行移動していくと、領域から出る直前に触れる点が最適です。平らな直線が凸多面体から離れる瞬間に触れているのは頂点（まれに辺全体）——だから線形計画の最適解は必ず頂点（端点）に存在します。 $n$ 変数なら、頂点は「 $n$ 本の制約が同時に等号で交わる点」。最適化とは「どの制約を等号で効かせるか（どの角か）」の組合せ選びに帰着します。

行列で書くと標準形は次の通りです。

\max_{x}\ c^\top x \quad \text{s.t.}\quad Ax\le b,\ \ x\ge 0

各資源制約に双対変数 $y_j\ge0$ を対応させた双対問題を考えると、

\min_{y}\ b^\top y \quad \text{s.t.}\quad A^\top y\ge c,\ \ y\ge 0

強双対定理により、最適では主問題の最大利益と双対問題の最小値が一致します。この最適な $y_j$ が資源 $j$ のシャドープライスで、定義は

y_j=\frac{\partial Z^\*}{\partial b_j}

——資源 $j$ の上限 $b_j$ を1単位ゆるめたときの最適利益の増分です。さらに補完スラック条件 $y_j\cdot(\text{slack}_j)=0$ が成り立ちます。すなわち、

資源が余っている（slack $_j>0$ ）なら $y_j=0$ ：余り物を足しても利益は増えない。
$y_j>0$ なら slack $_j=0$ ：その資源は**拘束（ボトルネック）**で、増やせば利益が伸びる。

これが「どの資源に投資すれば効くか」の答えになります。凸性・双対・KKT条件の証明と一般論は凸最適化の基礎へ（LP はそこで扱う凸最適化を線形に特殊化したもの）。本ノートはこの結果を linprog で計算し、増資源の再求解で裏取りします。

3. linprog で製品ミックスを解く（コード）

scipy.optimize.linprog は最小化ソルバです。利益最大化は目的係数の符号を反転（ $c\to-c$ ）して最小化に直し、得られた最小値の符号を戻します。2製品・3資源の例を解いて、最適生産量・最大利益・各資源の余り（slack）からボトルネックを判定し、シャドープライスを取り出します。

import numpy as np
from scipy.optimize import linprog

# 製品ミックス：製品A・製品B を何個つくると利益が最大か
# 利益（千円/個）：A=40, B=30
c_profit = np.array([40.0, 30.0])

# 資源制約 A_ub x <= b_ub
#   機械加工(時間/個)：2A + 1B <= 100
#   組立    (時間/個)：1A + 2B <=  80
#   原料    (kg/個)  ：1A + 1B <=  70
A_ub = np.array([
    [2.0, 1.0],   # 機械加工
    [1.0, 2.0],   # 組立
    [1.0, 1.0],   # 原料
])
b_ub = np.array([100.0, 80.0, 70.0])
res_names = ["機械加工(時間)", "組立(時間)", "原料(kg)"]

# linprog は「最小化」。利益最大化は目的を負にして最小化する（c -> -c）
res = linprog(c=-c_profit, A_ub=A_ub, b_ub=b_ub, bounds=[(0, None), (0, None)],
              method="highs")

x_opt = res.x
profit_opt = -res.fun                     # 符号を戻す
used = A_ub @ x_opt                        # 各資源の使用量
slack = b_ub - used                        # 余り（slack）
shadow = -res.ineqlin.marginals           # シャドープライス＝利益増/資源1単位

print("=== 最適生産計画 ===")
print(f"製品A の生産量 x_A = {x_opt[0]:.4f} 個")
print(f"製品B の生産量 x_B = {x_opt[1]:.4f} 個")
print(f"最大利益          = {profit_opt:.4f} 千円")
print()
print("=== 資源の使用状況（slack=0 が拘束制約＝ボトルネック）===")
print(f"{'資源':<14}{'使用量':>10}{'上限':>8}{'余りslack':>12}{'シャドープライス':>18}")
for i, name in enumerate(res_names):
    tag = "  <-拘束(ボトルネック)" if slack[i] < 1e-6 else "  （余裕あり）"
    print(f"{name:<14}{used[i]:>10.2f}{b_ub[i]:>8.1f}{slack[i]:>12.4f}{shadow[i]:>18.4f}{tag}")

出力：

=== 最適生産計画 ===
製品A の生産量 x_A = 40.0000 個
製品B の生産量 x_B = 20.0000 個
最大利益          = 2200.0000 千円

=== 資源の使用状況（slack=0 が拘束制約＝ボトルネック）===
資源                   使用量      上限     余りslack          シャドープライス
機械加工(時間)          100.00   100.0      0.0000           16.6667  <-拘束(ボトルネック)
組立(時間)             80.00    80.0      0.0000            6.6667  <-拘束(ボトルネック)
原料(kg)             60.00    70.0     10.0000            0.0000  （余裕あり）

出力の意味：最適は 製品A を40個・製品B を20個、最大利益 2200千円。注目は資源の使い方です——機械加工と組立は使用量が上限ぴったり（slack $=0$ ）で、ここが拘束制約＝ボトルネック。最適解はこの2本の制約が交わる頂点 $(40,20)$ に来ています（2変数なので2本が等号で効く）。一方原料は60kg使用で上限70kg、10kg余っています（slack $>0$ ）＝律速していません。シャドープライスを見ると、拘束している機械加工は 16.67千円/時間、組立は 6.67千円/時間 と正の値で、余っている原料は0。つまり「もし生産能力を増やすなら、原料を買い増しても1円も儲からない。効くのは機械加工の時間を増やすこと（1時間あたり16.67千円の利益増）」と読めます。ボトルネックがどこで、そこを1単位広げると利益がいくら伸びるか——LP はこれを同時に教えてくれます。

4. シャドープライスを再求解で検算し、実行可能領域を図示する（コード）

シャドープライスは「資源を1単位増やすと利益が $y_j$ 増える」という主張でした。本当かどうかは、上限 $b_j$ を実際に $+1$ して解き直し、利益の増分を見れば確かめられます。あわせて2製品なので、実行可能領域と最適頂点・等利益線を図示します。

import numpy as np
from scipy.optimize import linprog
import matplotlib.pyplot as plt
import japanize_matplotlib

c_profit = np.array([40.0, 30.0])
A_ub = np.array([[2.0, 1.0], [1.0, 2.0], [1.0, 1.0]])
b_ub = np.array([100.0, 80.0, 70.0])
res_names = ["機械加工", "組立", "原料"]

def solve(b):
    r = linprog(c=-c_profit, A_ub=A_ub, b_ub=b, bounds=[(0, None), (0, None)],
                method="highs")
    return -r.fun, r.x

base_profit, x_opt = solve(b_ub)
shadow = -linprog(c=-c_profit, A_ub=A_ub, b_ub=b_ub,
                  bounds=[(0, None), (0, None)], method="highs").ineqlin.marginals

print("=== シャドープライスの検算（資源を1単位増やして再求解）===")
print(f"基準の最大利益 = {base_profit:.4f} 千円")
print(f"{'資源':<8}{'シャドープライス':>16}{'+1単位後の利益':>16}{'実際の利益増':>14}")
for i, name in enumerate(res_names):
    b2 = b_ub.copy()
    b2[i] += 1.0
    p2, _ = solve(b2)
    print(f"{name:<8}{shadow[i]:>16.4f}{p2:>16.4f}{p2 - base_profit:>14.4f}")
print("-> 拘束資源では『利益増 = シャドープライス』、非拘束(原料)は増やしても利益不変")

# 図：実行可能領域（2製品）と最適頂点
xx = np.linspace(0, 60, 400)
plt.figure(figsize=(8.5, 7))
colors = ["#1f77b4", "#2ca02c", "#9467bd"]
for i, name in enumerate(res_names):
    a1, a2 = A_ub[i]
    yy = (b_ub[i] - a1 * xx) / a2          # 制約境界 a1*A + a2*B = b
    plt.plot(xx, yy, color=colors[i], lw=1.8,
             label=f"{name}: {a1:.0f}A+{a2:.0f}B<={b_ub[i]:.0f}")

# 実行可能領域（全制約を満たすグリッド点）を塗る
gx, gy = np.meshgrid(np.linspace(0, 60, 600), np.linspace(0, 60, 600))
feasible = ((A_ub[0, 0] * gx + A_ub[0, 1] * gy <= b_ub[0]) &
            (A_ub[1, 0] * gx + A_ub[1, 1] * gy <= b_ub[1]) &
            (A_ub[2, 0] * gx + A_ub[2, 1] * gy <= b_ub[2]))
plt.contourf(gx, gy, feasible.astype(float), levels=[0.5, 1.5],
             colors=["#ffd54f"], alpha=0.35)

# 等利益線（最適利益を通る）
P = base_profit
yiso = (P - c_profit[0] * xx) / c_profit[1]
plt.plot(xx, yiso, "--", color="#d62728", lw=1.6,
         label=f"等利益線 40A+30B={P:.0f}")

plt.plot(x_opt[0], x_opt[1], "*", color="black", ms=20,
         label=f"最適頂点 ({x_opt[0]:.0f}, {x_opt[1]:.0f})")
plt.xlabel("製品A の生産量"); plt.ylabel("製品B の生産量")
plt.title("線形計画：実行可能領域の頂点が最適（拘束2本の交点）")
plt.xlim(0, 60); plt.ylim(0, 60)
plt.legend(loc="upper right", fontsize=9)
plt.grid(alpha=0.3); plt.tight_layout()
plt.show()

出力：

=== シャドープライスの検算（資源を1単位増やして再求解）===
基準の最大利益 = 2200.0000 千円
資源              シャドープライス        +1単位後の利益        実際の利益増
機械加工             16.6667       2216.6667       16.6667
組立                6.6667       2206.6667        6.6667
原料                0.0000       2200.0000        0.0000
-> 拘束資源では『利益増 = シャドープライス』、非拘束(原料)は増やしても利益不変

出力の意味：機械加工を $+1$ 時間すると利益は $2200\to2216.67$ で $+16.67$ ——シャドープライス 16.67 とぴたり一致。組立 $+1$ は $+6.67$ でこれも一致。そして余っている原料を $+1$ kg しても利益は2200のまま（増分0）で、補完スラック条件（非拘束資源のシャドープライス $=0$ ）が数値で確かめられました。シャドープライスは机上の双対値ではなく、「ここを1単位広げたら本当にこれだけ儲かる」という運用上の限界価値なのです。図では、3本の制約境界が囲む黄色い実行可能領域の右上の角（機械加工＝青と組立＝緑の交点）で赤い等利益線が領域に接し、そこに最適頂点 $(40,20)$ の★が乗っています。原料（紫）はこの角より外側を通り、領域を縛っていない＝余っていることが視覚的にもわかります。利益を増やす向きへ平行移動した等利益線が、最後に触れるのは頂点——LP の幾何学的な正体です。

⚠️ よくある誤解

「生産量は連続値（分割可能）でよい」とは限らない：LP は $x_i$ が実数（小数）を取れる前提です。「製品Aを40.0個」のように割り切れれば問題ありませんが、「大型機械を2.7台」「フライトを3.4便」のように1個単位でしか作れないものは、 $x_i$ を整数に制約する整数計画（ILP）になります。ILP は LP よりはるかに難しく（一般にNP困難）、最適解は「LP解を丸めたもの」とは限りません。分割可能なら LP、不可分なら ILP と使い分けます。
「最適解は領域の真ん中あたりにある」ではない：線形目的の最適は必ず**頂点（角）**です。内部点が最適になることはありません（目的が線形ゆえ、必ずどこかの方向に動かせば改善できる）。だから LP は「どの制約を同時に効かせる角か」を選ぶ問題に帰着します。複数の頂点が同じ最適値を持つ（辺全体が最適）こともありますが、その場合も頂点に最適解が存在します。
「拘束していない資源を増やせば儲かる」ではない：余っている資源（slack $>0$ ）のシャドープライスは0。律速していないものをいくら買い増しても利益は1円も増えません。投資すべきは**拘束資源（ボトルネック）**だけ。ボトルネックを広げると、やがて律速が別の資源へ移り（ボトルネックとキャパシティの「ボトルネックは動く」）、そこからはシャドープライスが変わります。
「シャドープライスはどこまで増やしても一定」ではない：シャドープライスは局所的な値で、その資源を増やし続けると、いずれ別の制約が新たに拘束になり、最適頂点が切り替わってシャドープライスも変わります（有効範囲＝RHS感度分析）。「機械加工を1時間増やすと16.67千円」は今の頂点の近くでの話で、何十時間も増やせば原料が先に尽きてシャドープライスは下がります。一定なのは「ある範囲内」だけです。
「linprog にそのまま最大化を渡せる」ではない：scipy.optimize.linprog は最小化専用。利益最大化は目的係数を負にして渡し（ $c\to-c$ ）、得られた res.fun の符号を戻します。符号の戻し忘れは典型的なバグです。