21. 数理計画法
テーマ: 生産計画の線形計画問題と双対価格（シャドープライス）
学習点: scipy.optimize.linprog, 制約行列の組み立て, 双対問題の解釈, 感度分析
依存: NumPy, SciPy ／ 難易度: 上級
実行方法
uv run 21_linear_programming.py
スクリプト冒頭の PEP 723 メタデータ（# /// script）により、必要なライブラリは
uv が自動的に仮想環境へ導入します。事前の pip install は不要です。
解説
何をするプログラムか
工場で 3 つの製品 X・Y・Z を生産するとき、機械時間・労働時間・原材料・検査工程という 4 つの資源の保有量を超えない範囲で、限界利益の合計を最大にする生産量を決める――これは経営科学の古典である線形計画問題（LP）です。このスクリプトは SciPy の linprog で最適生産計画を求めたうえで、LP 特有の副産物であるシャドープライス（双対価格）を取り出します。
シャドープライスは「その資源をあと 1 単位調達できたら利益がいくら増えるか」を表す値で、資源の追加購入や残業の可否といった意思決定の判断材料になります。後半では双対問題を直接解いて強双対定理を確認し、右辺を動かす感度分析、さらに整数制約付きの求解まで行います。
コードの読みどころ
linprogは最小化専用なので、linprog(-profit, A_ub=A, b_ub=b, ...)と目的関数の符号を反転して渡し、結果は-res.funで最大利益に戻しています。- 制約行列
Aは「行 = 資源、列 = 製品」で組み立てられており、最適解の資源使用量はused = A @ res.xと行列積 1 行で計算できます。 - シャドープライスは
duals = -res.ineqlin.marginalsで取得します。HiGHS ソルバは符号を反転して返すため、マイナスを付けるのが約束です。 - 双対問題は
linprog(b, A_ub=-A.T, b_ub=-profit, ...)と、係数を転置・符号反転するだけで主問題と同じ関数で解けます。dres.funが主問題の最適値と一致することを出力で確認します。 - 感度分析のループでは
bb[0] += addと機械時間の右辺だけ増やして再求解し、利益の増分を 1 時間あたりに換算してシャドープライスと突き合わせます。 - 最後の
integrality=[1, 1, 1]を渡すだけで、同じlinprogが整数計画（MILP)として解いてくれます。
理論的背景
主問題 max c’x s.t. Ax≤b, x≥0 には、双対問題 min b’y s.t. A’y≥c, y≥0 が対応します。強双対定理により両者の最適値は一致し、双対変数 y_i は「制約 i の右辺 b_i を 1 単位緩めたときの目的関数の増加量」、すなわちシャドープライスです。また相補性条件から、余裕（スラック）が正の制約のシャドープライスは必ず 0 になります。使い切っていない資源を追加しても利益は増えない、という直観と一致します。
実行結果の見方
最適解は X=200・Y=400・Z=200 で総利益 3,080,000 円です。資源の表を見ると、機械時間・労働時間・原材料は余裕 0（ボトルネック）でシャドープライスが正、検査工程は余裕 700 が残るため 0 になっており、相補性条件どおりです。特に労働時間の 747.62 円が最も高く、増強すべき資源が一目で分かります。
感度分析では、機械時間を +50〜+400 時間まで増やす間は 1 時間あたりの利益増分が 390.48 円で一定（シャドープライスと一致）ですが、+800 時間では 385.71 円に下がります。基底が変わり別の制約が拘束的になるためで、シャドープライスが有効なのは一定範囲までという注意点が数値で確認できます。整数解は連続解と偶然一致し、整数ギャップは 0 円です。
ソースコード
# /// script
# requires-python = ">=3.11"
# dependencies = [
# "numpy",
# "scipy",
# ]
# ///
"""21: 数理計画法 ---------------------------------------------------------
テーマ: 生産計画の線形計画問題と双対価格（シャドープライス）
学習点: scipy.optimize.linprog, 制約行列の組み立て, 双対問題の解釈,
感度分析
根拠: 主問題 max c'x s.t. Ax<=b, x>=0 に対し双対問題 min b'y s.t. A'y>=c, y>=0。
強双対定理より最適値は一致し、双対変数 y_i は制約 i の右辺を1単位
緩めたときの目的関数の増加量（シャドープライス）を表す。
相補性条件より、余裕のある制約のシャドープライスは 0 になる。
"""
import numpy as np
from scipy.optimize import linprog
def main() -> None:
products = ["製品X", "製品Y", "製品Z"]
profit = np.array([4200.0, 3100.0, 5000.0]) # 単位あたり限界利益
resources = ["機械時間", "労働時間", "原材料kg", "検査工程"]
# 各製品が消費する資源量（行=資源, 列=製品）
A = np.array([
[3.0, 2.0, 5.0],
[4.0, 3.0, 4.0],
[2.0, 4.0, 3.0],
[1.0, 1.0, 1.0],
])
# 検査工程は十分な余裕がある（＝ボトルネックでない）ように設定
b = np.array([2400.0, 2800.0, 2600.0, 1500.0]) # 利用可能量
# linprog は最小化なので目的関数の符号を反転する
res = linprog(-profit, A_ub=A, b_ub=b, bounds=[(0, None)] * 3,
method="highs")
assert res.success, res.message
print("=== 最適生産計画 ===")
for p, q in zip(products, res.x):
print(f" {p}: {q:>9.3f} 単位 (限界利益 {profit[list(products).index(p)]:,.0f}円)")
print(f" 最大総利益 = {-res.fun:,.2f} 円\n")
print(f"{'資源':<10}{'使用量':>11}{'保有量':>11}{'余裕':>11}"
f"{'シャドープライス':>18}")
print("-" * 62)
used = A @ res.x
duals = -res.ineqlin.marginals # highs は符号が反転して返る
for r, u, cap, y in zip(resources, used, b, duals):
print(f"{r:<10}{u:>11.2f}{cap:>11.1f}{cap - u:>11.3f}{y:>18,.2f}")
print("\n 解釈: シャドープライスは、その資源を1単位追加調達したときに")
print(" 増える利益の上限。これを超える価格で買うのは不合理。")
print(" 余裕(スラック)が正の制約はボトルネックではないため価格0。\n")
print("[双対問題を直接解いて強双対定理を確認]")
dres = linprog(b, A_ub=-A.T, b_ub=-profit,
bounds=[(0, None)] * len(b), method="highs")
print(f" 主問題の最適値 = {-res.fun:,.4f}")
print(f" 双対問題の最適値 = {dres.fun:,.4f} -> 一致（強双対定理）")
print(f" 双対解 y = {np.round(dres.x, 4)}\n")
print("[右辺の感度分析] 機械時間を増やしたときの利益")
base = -res.fun
for add in [0, 50, 100, 200, 400, 800]:
bb = b.copy()
bb[0] += add
r2 = linprog(-profit, A_ub=A, b_ub=bb, bounds=[(0, None)] * 3,
method="highs")
inc = -r2.fun - base
print(f" +{add:>4.0f}時間: 利益 {-r2.fun:>12,.1f}円 "
f"(増分 {inc:>10,.1f} / 1時間あたり {inc/add if add else 0:>8,.2f})")
print(" ※ シャドープライスが一定なのは基底が変わらない範囲まで。"
"\n 増やしすぎると別の制約が拘束的になり限界価値が下がる（区分線形）。\n")
print("[整数制約がある場合] 製品は整数単位でしか作れないなら")
ires = linprog(-profit, A_ub=A, b_ub=b, bounds=[(0, None)] * 3,
integrality=[1, 1, 1], method="highs")
print(f" 整数解: {np.round(ires.x, 0)} -> 利益 {-ires.fun:,.0f}円")
print(f" 連続緩和との差（整数ギャップ）= {base + ires.fun:,.2f}円")
if __name__ == "__main__":
main()
実行結果
=== 最適生産計画 ===
製品X: 200.000 単位 (限界利益 4,200円)
製品Y: 400.000 単位 (限界利益 3,100円)
製品Z: 200.000 単位 (限界利益 5,000円)
最大総利益 = 3,080,000.00 円
資源 使用量 保有量 余裕 シャドープライス
--------------------------------------------------------------
機械時間 2400.00 2400.0 0.000 390.48
労働時間 2800.00 2800.0 0.000 747.62
原材料kg 2600.00 2600.0 0.000 19.05
検査工程 800.00 1500.0 700.000 0.00
解釈: シャドープライスは、その資源を1単位追加調達したときに
増える利益の上限。これを超える価格で買うのは不合理。
余裕(スラック)が正の制約はボトルネックではないため価格0。
[双対問題を直接解いて強双対定理を確認]
主問題の最適値 = 3,080,000.0000
双対問題の最適値 = 3,080,000.0000 -> 一致（強双対定理）
双対解 y = [390.4762 747.619 19.0476 0. ]
[右辺の感度分析] 機械時間を増やしたときの利益
+ 0時間: 利益 3,080,000.0円 (増分 0.0 / 1時間あたり 0.00)
+ 50時間: 利益 3,099,523.8円 (増分 19,523.8 / 1時間あたり 390.48)
+ 100時間: 利益 3,119,047.6円 (増分 39,047.6 / 1時間あたり 390.48)
+ 200時間: 利益 3,158,095.2円 (増分 78,095.2 / 1時間あたり 390.48)
+ 400時間: 利益 3,236,190.5円 (増分 156,190.5 / 1時間あたり 390.48)
+ 800時間: 利益 3,388,571.4円 (増分 308,571.4 / 1時間あたり 385.71)
※ シャドープライスが一定なのは基底が変わらない範囲まで。
増やしすぎると別の制約が拘束的になり限界価値が下がる（区分線形）。
[整数制約がある場合] 製品は整数単位でしか作れないなら
整数解: [200. 400. 200.] -> 利益 3,080,000円
連続緩和との差（整数ギャップ）= 0.00円