#!/usr/bin/env python3 """Generate all figures for the Cyber Blast Radius article from results/*.csv|json. Run after cbr_model.py, grid_sim.py, correlated_failure.py, monte_carlo.py: python3 model/make_figures.py → results/figures/fig*.png + captions.md """ from __future__ import annotations import csv import json import os import tempfile from pathlib import Path os.environ.setdefault("MPLCONFIGDIR", tempfile.mkdtemp(prefix="mpl-")) import matplotlib # noqa: E402 matplotlib.use("Agg") import matplotlib.pyplot as plt # noqa: E402 import numpy as np # noqa: E402 from matplotlib.patches import FancyArrowPatch, FancyBboxPatch # noqa: E402 ROOT = Path(__file__).resolve().parents[1] RES = ROOT / "results" FIG = RES / "figures" FIG.mkdir(parents=True, exist_ok=True) plt.rcParams["font.family"] = ["Hiragino Sans", "Noto Sans CJK JP", "sans-serif"] plt.rcParams["axes.unicode_minus"] = False plt.rcParams["figure.dpi"] = 150 plt.rcParams["savefig.dpi"] = 170 plt.rcParams["axes.grid"] = True plt.rcParams["grid.alpha"] = 0.3 C_A = "#b23a48" # centralized / red C_B = "#1f6f8b" # distributed / blue C_H = "#d98e04" # hierarchical / amber C_F = "#5c8d5a" # federated / green C_G = "#666666" captions: list[tuple[str, str]] = [] def rows(name: str) -> list[dict]: with (RES / name).open(encoding="utf-8") as f: return list(csv.DictReader(f)) def fnum(v) -> float: try: return float(v) except (TypeError, ValueError): return float("nan") def save(fig, name: str, caption: str) -> None: fig.tight_layout() fig.savefig(FIG / name, bbox_inches="tight", facecolor="white") plt.close(fig) captions.append((name, caption)) print("wrote", name) def box(ax, x, y, w, h, text, fc, fs=8, ec="#333", lw=1.0): ax.add_patch(FancyBboxPatch((x, y), w, h, boxstyle="round,pad=0.02,rounding_size=0.03", fc=fc, ec=ec, lw=lw)) ax.text(x + w / 2, y + h / 2, text, ha="center", va="center", fontsize=fs, wrap=True) def arrow(ax, x1, y1, x2, y2, color="#333", style="-|>", lw=1.0, ls="-"): ax.add_patch(FancyArrowPatch((x1, y1), (x2, y2), arrowstyle=style, color=color, lw=lw, linestyle=ls, mutation_scale=10)) # ---------------------------------------------------------------- fig01 architectures def fig01_architectures(): fig, axes = plt.subplots(1, 3, figsize=(13, 5.2)) for ax in axes: ax.set_xlim(0, 10); ax.set_ylim(0, 10); ax.axis("off"); ax.grid(False) # A: centralized ax = axes[0] ax.set_title("モデル A：集中制御（単一制御面）", fontsize=10, color=C_A) box(ax, 2.5, 8.2, 5, 1.3, "全国DERMS / クラウド\n（1つの認証情報・1つのAPI）", "#f7d9dc", fs=8) for i, x in enumerate(np.linspace(0.8, 8.4, 6)): box(ax, x, 3.6, 1.0, 1.0, "PCS", "#fbeaec", fs=7) arrow(ax, 5, 8.2, x + 0.5, 4.6, color=C_A) box(ax, 0.6, 1.2, 8.8, 1.3, "命令はそのまま実行（ローカル検証なし）\nCBR = 到達可能な全容量", "#fff", fs=8, ec=C_A) ax.text(5, 0.4, "侵害1件 → 35 GW が同時に動く", ha="center", fontsize=8, color=C_A) # H: hierarchical / federated ax = axes[1] ax.set_title("モデル H/F：階層・連邦（独立ドメイン）", fontsize=10, color=C_F) box(ax, 3.2, 8.4, 3.6, 1.1, "系統運用者（助言・市場）", "#e7efe6", fs=8) for j, x in enumerate((0.6, 3.7, 6.8)): box(ax, x, 6.0, 2.6, 1.1, f"ドメイン {j+1}\n独立認証・権限上限", "#dfe9de", fs=7) arrow(ax, 5, 8.4, x + 1.3, 7.1, color=C_F, ls="--") for k in range(3): xx = x + 0.1 + k * 0.85 box(ax, xx, 3.8, 0.7, 0.9, "PCS", "#eef4ee", fs=6) arrow(ax, x + 1.3, 6.0, xx + 0.35, 4.7, color=C_F) box(ax, 0.6, 1.2, 8.8, 1.3, "侵害はドメイン内に閉じる（認証が独立なら）\nCBR = ドメイン容量 ≤ 権限上限", "#fff", fs=8, ec=C_F) ax.text(5, 0.4, "100 ドメイン → 侵害1件で 350 MW", ha="center", fontsize=8, color=C_F) # B: distributed autonomous ax = axes[2] ax.set_title("モデル B：分散自律（ローカル制約エンジン）", fontsize=10, color=C_B) box(ax, 3.0, 8.4, 4.0, 1.1, "上位通信＝助言・最適化層\n（任意・切れても動く）", "#d9e8ee", fs=7) for i, x in enumerate((0.4, 3.6, 6.8)): box(ax, x, 5.3, 2.8, 2.2, "PCS\n局所観測\nV / f / RoCoF / SOC\n制約エンジン\n受理・修正・拒否", "#e6f0f4", fs=6.2) arrow(ax, 5, 8.4, x + 1.4, 7.5, color=C_B, ls="--") box(ax, 0.6, 2.9, 8.8, 1.6, "コマンド・エンベロープ：ΔP/Δt・ΔPmax・N_max・P_region,max\n上位命令はこの枠を超えない（ρ：ローカル強制の生存率）", "#fff", fs=7.5, ec=C_B) box(ax, 0.6, 1.0, 8.8, 1.3, "物理（周波数・電圧）がソフトの権限に優先\nCBR_eff = reach × min(CBR, A_max) × (1 − ρ(1 − e))", "#fff", fs=7.5, ec=C_B) ax.text(5, 0.3, "100 ドメイン＋エンベロープ → 117 MW", ha="center", fontsize=8, color=C_B) save(fig, "fig01_architectures.png", "図1　三つのアーキテクチャ。左：単一制御面（命令はそのまま実行）。中：独立ドメインに分割（侵害はドメイン内に閉じる）。右：分散自律（各PCSが局所観測で上位命令を受理・修正・拒否し、エンベロープを超える変化を物理的に拒む）。") # ---------------------------------------------------------------- fig02 CBR decomposition def fig02_cbr_decomposition(): r = rows("cbr_decomposition.csv") labels = [x["unit"].replace("CBR_", "") for x in r] raw = [fnum(x["cbr_raw_mw"]) for x in r] noauto = [fnum(x["cbr_eff_no_autonomy_mw"]) for x in r] auto = [fnum(x["cbr_eff_with_autonomy_mw"]) for x in r] fig, ax = plt.subplots(figsize=(8.5, 4.6)) xx = np.arange(len(labels)); w = 0.27 ax.bar(xx - w, raw, w, label="CBR_raw（設置容量）", color="#bbb") ax.bar(xx, noauto, w, label="CBR_eff 自律なし（reach=0.7）", color=C_A) ax.bar(xx + w, auto, w, label="CBR_eff 自律あり（エンベロープ30%・ρ=0.95）", color=C_B) ax.set_yscale("log"); ax.set_ylabel("MW（対数）") ax.set_xticks(xx); ax.set_xticklabels(labels) fcr = fnum(rows("architecture_comparison.csv")[0]["fcr_mw"]) ax.axhline(fcr, color="k", ls="--", lw=0.8); ax.text(len(labels) - 0.5, fcr * 1.15, f"一次調整力 FCR ≈ {fcr:,.0f} MW", ha="right", fontsize=8) for i in range(len(labels)): ax.annotate(f"{auto[i]:,.1f}" if auto[i] < 100 else f"{auto[i]:,.0f}", (xx[i] + w, auto[i]), ha="center", va="bottom", fontsize=7, color=C_B) ax.annotate(f"{noauto[i]:,.0f}" if noauto[i] >= 1 else f"{noauto[i]:.4f}", (xx[i], noauto[i]), ha="center", va="bottom", fontsize=7, color=C_A) ax.legend(fontsize=8, loc="upper left") save(fig, "fig02_cbr_decomposition.png", "図2　Cyber Blast Radius の分解（基準ケース：1,000万台×5 kW＝50 GW）。侵害単位を機器→アグリゲータ→ベンダー→クラウド→地域と広げたときの到達容量。破線は系統の一次調整力。クラウド1件の侵害は自律なしで 35 GW、エンベロープ付きでも 11.7 GW と調整力を大きく超える。") # ---------------------------------------------------------------- fig03 severity calibration def fig03_severity(): cal = rows("severity_calibration.csv") x = np.linspace(0, 6, 300) fig, ax = plt.subplots(figsize=(8.5, 4.6)) for k, x0, ls, lab in ((3, 1.7, "-", "基準 k=3, x₀=1.7"), (2, 1.7, "--", "k=2"), (4, 1.7, ":", "k=4"), (3, 1.3, "-.", "x₀=1.3（脆弱な系統）")): ax.plot(x, 1 - np.exp(-(x / x0) ** k), ls=ls, color=C_G if k != 3 or x0 != 1.7 else "k", lw=1.6 if (k == 3 and x0 == 1.7) else 1.0, label=lab) for c in cal: ax.scatter(fnum(c["x"]), fnum(c["expected_order"]), s=40, color=C_A, zorder=5) ax.annotate(c["case"], (fnum(c["x"]), fnum(c["expected_order"])), textcoords="offset points", xytext=(6, -10 if "Hokkaido" in c["case"] else 6), fontsize=7) sw = rows("grid_step_sweep.csv") xs = [fnum(s["x_over_fcr"]) for s in sw]; ufls = [1.0 if "UFLS" in s["class"] else (0.5 if "degraded" in s["class"] else 0.0) for s in sw] ax.step(xs, ufls, where="post", color=C_B, lw=1, alpha=0.6, label="スイング方程式モデル（基準系統）：0=収束 0.5=劣化 1=UFLS") ax.set_xlabel("x = 同時擾乱 ÷ 系統の一次調整力（FCR）"); ax.set_ylabel("重大度 S(x) ＝ P(広域障害 | 擾乱)") ax.set_ylim(-0.03, 1.05); ax.legend(fontsize=7.5, loc="lower right") save(fig, "fig03_severity_calibration.png", "図3　重大度関数の較正。横軸は擾乱を一次調整力で割った無次元量 x。赤点は実事故（仏2023年 2.66 GW 無UFLS、ENTSO-E基準事故 3 GW、北海道2018 など）に対する「この規模が広域障害になった／ならなかった」の経験的目安、黒線はワイブル型の採用関数、青階段は本稿のスイング方程式モデルの判定。x≈2 を超えると UFLS が始まる。") # ---------------------------------------------------------------- fig04 architecture comparison def fig04_architecture_comparison(): r = rows("architecture_comparison.csv") order = ["A0_centralized_baseline", "A_centralized_certified", "A_centralized_certified_10x", "A+_centralized_server_cap", "H_hierarchical_10", "F_federated_100", "C_cellular_10000", "B4_single_plane_envelope_no_cap", "B1_distributed_envelope_no_cap", "B_distributed_autonomous", "B3_distributed_certified", "B5_distributed_hw_envelope", "B6_distributed_hw_vendor_cap12", "B7_distributed_hw_vendor_cap5"] short = {"A0_centralized_baseline": "A0 集中", "A_centralized_certified": "A 集中＋認証(p/3)", "A_centralized_certified_10x": "A 集中＋認証(p/10)", "A+_centralized_server_cap": "A+ 集中＋ｻｰﾊﾞｰ側上限", "H_hierarchical_10": "H 階層10", "F_federated_100": "F 連邦100", "C_cellular_10000": "C セル10,000", "B4_single_plane_envelope_no_cap": "B4 単一面＋ｴﾝﾍﾞﾛｰﾌﾟ", "B1_distributed_envelope_no_cap": "B1 100ﾄﾞﾒｲﾝ＋ｴﾝﾍﾞﾛｰﾌﾟ", "B_distributed_autonomous": "B 分散自律(ｿﾌﾄのみ)", "B3_distributed_certified": "B＋認証(p/3)", "B5_distributed_hw_envelope": "B5 ＋HW強制ｴﾝﾍﾞﾛｰﾌﾟ", "B6_distributed_hw_vendor_cap12": "B6 ＋ﾍﾞﾝﾀﾞｰ上限12%", "B7_distributed_hw_vendor_cap5": "B7 ＋ﾍﾞﾝﾀﾞｰ上限5%"} d = {x["architecture"]: x for x in r} names = [short[k] for k in order] cbr = [fnum(d[k]["cbr_eff_domain_mw"]) for k in order] cbr_v = [fnum(d[k]["cbr_eff_vendor_mw"]) for k in order] plane = [fnum(d[k]["fleet_systemic_risk"]) for k in order] vend = [fnum(d[k]["fleet_systemic_risk_vendor_channel"]) for k in order] total = [fnum(d[k]["fleet_systemic_risk_total"]) for k in order] colors = [C_A if k.startswith("A") else C_H if k.startswith("H") else C_F if k.startswith(("F", "C")) else C_B for k in order] fig, axes = plt.subplots(1, 3, figsize=(14, 5.6)) y = np.arange(len(names)); h = 0.38 ax = axes[0] ax.barh(y - h / 2, cbr, h, color=colors, label="制御面／ドメイン1件の侵害") ax.barh(y + h / 2, cbr_v, h, color=colors, alpha=0.45, hatch="//", label="首位ベンダー1社の侵害（クラウド／FW）") ax.set_yticks(y); ax.set_yticklabels(names, fontsize=8); ax.set_xscale("log"); ax.invert_yaxis(); ax.set_title("最大 Cyber Blast Radius（MW）", fontsize=9) fcr = fnum(r[0]["fcr_mw"]); ax.axvline(fcr, color="k", ls="--", lw=0.8); ax.text(fcr, -0.9, "FCR", fontsize=7, ha="center"); ax.legend(fontsize=7, loc="lower right") ax = axes[1] ax.barh(y - h / 2, plane, h, color=colors, label="制御面経路（ドメイン数ぶんの攻撃面）") ax.barh(y + h / 2, vend, h, color=colors, alpha=0.45, hatch="//", label="ベンダー経路（全ベンダー）") ax.set_yticks(y); ax.set_yticklabels([]); ax.set_xscale("log"); ax.invert_yaxis(); ax.set_title("系統リスク 経路別（広域障害・件/年、フリート全体）", fontsize=9); ax.legend(fontsize=7, loc="lower left"); ax.set_xlim(1e-7, 0.5) ax = axes[2] ax.barh(y, total, 0.6, color=colors) for i, v in enumerate(total): ax.annotate(f"{v:.1e}", (v, i), textcoords="offset points", xytext=(4, 0), va="center", fontsize=7) ax.set_yticks(y); ax.set_yticklabels([]); ax.set_xscale("log"); ax.invert_yaxis(); ax.set_title("系統リスク 合計（両経路）", fontsize=9); ax.set_xlim(1e-4, 0.5) ax.axvline(total[1], color=C_A, ls=":", lw=0.9); ax.text(total[1], len(names) - 0.4, "A＋認証", fontsize=7, color=C_A, ha="center") save(fig, "fig04_architecture_comparison.png", "図4　14 構成の比較（フリート全体・年間）。左：侵害1件で同時に動く最大容量（実線＝制御面／1ドメイン、斜線＝首位ベンダーのファームウェア系統）。中：フリート全体の年間広域障害確率（UFLS 以上）を、制御面経路（n ドメイン分の攻撃面と共通原因 β を含む）とベンダー経路（クラウド／OTA／ファームウェア）に分けたもの。右：合計。ドメイン分割とエンベロープは制御面経路を 1〜2 桁下げるが、ベンダー経路は権限上限を迂回し、ソフトのみのエンベロープはファームウェアと一緒に消えるため、B（ソフトのみ）の合計は A＋認証と同程度にとどまる。合計を下げるのは独立監視による出力下限の保持（B5）とファームウェア系統の到達上限（B6・B7）であり、集中側でもサーバー側の集計上限（A+）が同じ桁まで下げる。") # ---------------------------------------------------------------- fig05 attack traces def fig05_attack_traces(): tr = rows("attack_traces.csv") sc = rows("attack_scenarios.csv") by = {} for t in tr: by.setdefault(t["scenario"], []).append((fnum(t["t_s"]), fnum(t["f_hz"]), fnum(t["attack_mw"]))) fig, axes = plt.subplots(1, 2, figsize=(13, 4.8)) pick = [("A1 cloud compromise, national plane, 50% cut", C_A, "-", "A1 集中：50%削減 17.5 GW"), ("A3 vendor cloud compromise (largest vendor), full OFF", C_A, "--", "A3 最大ベンダークラウド 10.5 GW"), ("H1 one regional domain of 10, full OFF", C_H, "-", "H1 10ドメイン中1つ 3.5 GW"), ("B2 national plane compromised, autonomy ON (no cap)", C_B, ":", "B2 単一面＋エンベロープ 11.7 GW"), ("B5b largest vendor firmware compromised, independent-monitor floor clamp survives (ramp lost)", C_B, "--", "B5b ベンダーFW侵害・下限クランプ残存 3.5 GW（ランプ消失）"), ("B5c firmware-line reach capped at 12% of fleet, floor clamp survives", C_B, "-.", "B5c FW系統の到達 12% 上限 1.4 GW"), ("A4 single plane with server-side cap 1,500 MW, API credential compromise", C_A, "-.", "A4 集中＋サーバー側上限 1.5 GW"), ("F1 one domain of 100, full OFF", C_F, "-", "F1 100ドメイン中1つ 350 MW"), ("B1 one domain of 100, autonomy ON", C_B, "-", "B1 分散自律 117 MW")] ax = axes[0] for key, col, ls, lab in pick: pts = by.get(key, []) if not pts: continue t = np.array([p[0] for p in pts]); f = np.array([p[1] for p in pts]) ax.plot(t, f, color=col, ls=ls, lw=1.4, label=lab) ax.axhline(48.5, color="k", lw=0.6, ls=":"); ax.text(1, 48.52, "UFLS 第1段 48.5 Hz（基準系統の仮定）", fontsize=7) ax.axhline(47.5, color="k", lw=0.6, ls=":"); ax.text(1, 47.52, "発電機脱落 47.5 Hz", fontsize=7) ax.set_xlim(0, 60); ax.set_ylim(47.3, 50.2); ax.set_xlabel("秒"); ax.set_ylabel("系統周波数 Hz"); ax.set_title("同一攻撃（クラウド侵害）に対する周波数応答：最初の60秒", fontsize=9) ax.legend(fontsize=7, loc="lower right") ax = axes[1] labels = [s["scenario"].split(" ")[0] for s in sc]; aff = [fnum(s["affected_mw"]) for s in sc]; nad = [fnum(s["nadir_hz"]) for s in sc] cols = [C_A if l.startswith("A") else C_H if l.startswith("H") else C_F if l.startswith("F") else C_B for l in labels] ax.scatter(aff, nad, c=cols, s=50, zorder=5) for l, a, n_, s in zip(labels, aff, nad, sc): ax.annotate(f"{l} ({s['class'].split(' ')[0]})", (a, n_), textcoords="offset points", xytext=(5, 4), fontsize=7) ax.set_xscale("log"); ax.set_xlabel("同時に失われた DER 出力 MW（対数）"); ax.set_ylabel("周波数最低値 Hz"); ax.set_title("影響容量 vs 周波数最低値（全シナリオ）", fontsize=9) ax.axhline(48.5, color="k", lw=0.6, ls=":"); ax.axvline(1350, color="k", lw=0.6, ls="--"); ax.text(1400, 50.05, "FCR", fontsize=7); ax.axvline(3600, color="k", lw=0.6, ls="-."); ax.text(3700, 50.05, "FCR+FRR", fontsize=7) save(fig, "fig05_attack_traces.png", "図5　同一の攻撃（制御クラウドまたはファームウェア系統の侵害で到達可能な DER を一斉削減）を単一エリアのスイング方程式モデルに与えたときの周波数応答。左：代表シナリオの最初の60秒。右：全シナリオの影響容量と周波数最低値。A1・A3・H1・B2・B5a・B5b は負荷遮断（UFLS）に至り、F1・B1・B4 は調整力の範囲に収まる。B5b では出力下限（30%）は残るがランプ制限はファームウェアと一緒に失われるため 3.5 GW がステップで落ちる。ファームウェア系統の到達を 12% に抑えた B5c と、サーバー側に 1,500 MW の集計上限を置いた A4 は、どちらも調整力をわずかに超える「劣化」域に収まる。") # ---------------------------------------------------------------- fig06 break-even heatmap def fig06_break_even(): g = rows("break_even_grid.csv") pf = sorted({fnum(r["p_reduction_factor"]) for r in g}); cm = sorted({fnum(r["blast_radius_multiplier"]) for r in g}) Z = np.full((len(pf), len(cm)), np.nan); L = np.full((len(pf), len(cm)), np.nan) for r in g: i = pf.index(fnum(r["p_reduction_factor"])); j = cm.index(fnum(r["blast_radius_multiplier"])) Z[i, j] = np.log10(max(fnum(r["risk_ratio_A_over_B"]), 1e-6)); L[i, j] = 1.0 if r["centralized_safer_linear"] == "True" else 0.0 fig, axes = plt.subplots(1, 2, figsize=(12.5, 4.8)) ax = axes[0] im = ax.imshow(Z, origin="lower", cmap="RdBu_r", vmin=-4, vmax=4, aspect="auto") ax.set_xticks(range(len(cm))); ax.set_xticklabels([f"×{c:g}" for c in cm], fontsize=7); ax.set_yticks(range(len(pf))); ax.set_yticklabels([f"1/{p:g}" for p in pf], fontsize=7) ax.set_xlabel("集中化による Blast Radius の倍率（基準 10 MW に対して）"); ax.set_ylabel("認証・集中管理による侵害確率の低減") ax.set_title("log10（集中リスク ÷ 分散リスク）：赤＝集中の方が危険", fontsize=9); ax.grid(False) for i in range(len(pf)): for j in range(len(cm)): ax.text(j, i, f"{Z[i,j]:+.1f}", ha="center", va="center", fontsize=6, color="w" if abs(Z[i, j]) > 2 else "k") fig.colorbar(im, ax=ax, shrink=0.8) ax = axes[1] k = rows("required_p_reduction_by_k.csv") for kk, col in ((1, "#999"), (2, C_F), (3, "k"), (4, C_B)): pts = [(fnum(r["cbr_a_mw"]), fnum(r["p_reduction_needed"])) for r in k if fnum(r["severity_k"]) == kk] pts.sort(); ax.plot([p[0] for p in pts], [min(p[1], 1e8) for p in pts], marker="o", ms=3, color=col, label=f"k={kk}" + ("（線形）" if kk == 1 else "（基準）" if kk == 3 else "")) ax.set_xscale("log"); ax.set_yscale("log"); ax.set_xlabel("集中側の Blast Radius MW（分散側 10 MW）"); ax.set_ylabel("リスクを等しくするのに必要な侵害確率の低減倍率") ax.axvline(1350, color="k", ls="--", lw=0.8); ax.text(1400, 2, "FCR", fontsize=7); ax.set_title("損益分岐：Blast Radius を広げた分を確率で埋めるには", fontsize=9); ax.legend(fontsize=8) save(fig, "fig06_break_even.png", "図6　損益分岐分析。左：侵害確率を 1/p に下げる代わりに Blast Radius を ×Y に広げたときのリスク比（対数）。確率を 1/10 にしても Blast Radius が ×100 なら集中側が危険。右：集中側の Blast Radius に対して、分散側（10 MW）とリスクを等しくするために必要な確率低減倍率。重大度関数の指数 k が結論を左右するため k=1〜4 を併記した。調整力（破線）を超えると、確率をどれだけ下げても重大度は飽和し、比較は確率だけの勝負になる。") # ---------------------------------------------------------------- fig07 CCF def fig07_ccf(): fl = rows("ccf_fleets.csv"); summ = json.loads((RES / "ccf_summary.json").read_text()) fig, axes = plt.subplots(1, 2, figsize=(13, 4.8)) ax = axes[0] short = {"monoculture (1 vendor, 1 cloud, 1 firmware, no autonomy)": "単一ベンダー・自律なし", "diverse (8 vendors, 100 domains, no autonomy)": "8ベンダー・100ドメイン・自律なし", "diverse + software envelope (lost on firmware compromise)": "同＋ソフトｴﾝﾍﾞﾛｰﾌﾟ", "diverse + hardware-enforced envelope": "同＋HW強制ｴﾝﾍﾞﾛｰﾌﾟ", "monoculture + hardware-enforced envelope": "単一ベンダー＋HW強制", "very diverse (20 vendors, top 15%) + hardware envelope": "20ベンダー(首位15%)＋HW強制"} cols = [C_A, C_H, C_B, C_B, C_A, C_F] for (name, h), col in zip(summ["histograms"].items(), cols): e = np.array(h["edges_mw"]); c = np.array(h["counts"], dtype=float); c = c / c.sum() tail = np.cumsum(c[::-1])[::-1] ax.step(e[:-1], tail, where="post", color=col, lw=1.3, label=short.get(name, name), ls="-" if "hardware" not in name else "--") ax.axvline(summ["fcr_mw"], color="k", ls="--", lw=0.8); ax.text(summ["fcr_mw"] * 1.1, 0.5, "FCR", fontsize=7) ax.set_xscale("log"); ax.set_yscale("log"); ax.set_ylim(1e-4, 1.2); ax.set_xlim(1, 6e4) ax.set_xlabel("1年に同時に敵の制御下に入る容量 MW"); ax.set_ylabel("超過確率（年あたり）"); ax.set_title("共通原因故障モンテカルロ：年間影響容量の裾", fontsize=9); ax.legend(fontsize=6.5, loc="lower left") ax = axes[1] sw = rows("ccf_vendor_sweep.csv") for hw, col, lab in (("False", C_A, "自律なし"), ("True", C_B, "HW強制エンベロープ")): pts = [(fnum(r["n_vendors"]), fnum(r["p_exceed_fcr_per_year"]), fnum(r["p99_mw"])) for r in sw if r["hardware_envelope"] == hw] ax.plot([p[0] for p in pts], [p[1] for p in pts], marker="o", color=col, label=f"P(影響>FCR)/年 — {lab}") ax.set_xlabel("ベンダー数（首位シェアは 1/n〜30%）"); ax.set_ylabel("年に一度でも調整力を超える確率"); ax.set_title("多様性のパラドクス：ベンダーを増やすと「どこかが侵害される」確率は上がる", fontsize=9) ax.legend(fontsize=7.5) ax2 = ax.twinx(); ax2.grid(False) for hw, col in (("False", C_A), ("True", C_B)): pts = [(fnum(r["n_vendors"]), fnum(r["p99_mw"])) for r in sw if r["hardware_envelope"] == hw] ax2.plot([p[0] for p in pts], [p[1] for p in pts], color=col, ls=":", marker="s", ms=3) ax2.set_yscale("log"); ax2.set_ylabel("99パーセンタイル影響容量 MW（点線）") save(fig, "fig07_correlated_failure.png", "図7　共通原因故障（ベンダークラウド・ファームウェア・アグリゲータ認証を共有する群が一斉に落ちる）のモンテカルロ。左：年間の同時影響容量の超過確率。単一ベンダー（赤実線）は確率こそ低いが裾が 35 GW まで伸びる。右：ベンダー数を増やすと「どこかのベンダーが侵害される」確率は上がる一方、99パーセンタイル容量は下がる。多様性は「侵害の頻度」でなく「一回あたりの大きさ」を抑える対策であり、各ベンダーのシェア上限とエンベロープを伴って初めて調整力の範囲に収まる。") # ---------------------------------------------------------------- fig08 Monte Carlo def fig08_monte_carlo(): mc = rows("mc_samples.csv"); summ = json.loads((RES / "monte_carlo_summary.json").read_text()) ra = np.array([fnum(r["risk_A"]) for r in mc]); rb = np.array([fnum(r["risk_B"]) for r in mc]); rf = np.array([fnum(r["risk_Bfull"]) for r in mc]) rp = np.array([fnum(r["risk_Aplus"]) for r in mc]) pa = np.array([fnum(r["plane_A"]) for r in mc]); pb = np.array([fnum(r["plane_B"]) for r in mc]) fig, axes = plt.subplots(1, 4, figsize=(17, 4.6)) lo, hi = 1e-9, 1 T = summ["total"]; P = summ["plane_channel_only"] for ax, xa, yb, title, col, xl in ( (axes[0], pa, pb, f"制御面経路のみ：B が低い {P['B_vs_A']['share_lower']:.0%}、1/10 以下 {P['B_vs_A']['share_lower_by_10x']:.0%}", C_B, "集中 A の制御面リスク /年"), (axes[1], ra, rb, f"合計：B（ソフトのみ）が A より低い {T['B_vs_A']['share_lower']:.0%}、中央値比 {T['B_vs_A']['median_ratio']:.1f}×", C_G, "集中 A の合計リスク /年"), (axes[2], ra, rf, f"合計：B＋独立監視＋FW到達12% が A より低い {T['Bfull_vs_A']['share_lower']:.0%}、中央値比 {T['Bfull_vs_A']['median_ratio']:.0f}×", C_F, "集中 A の合計リスク /年"), (axes[3], rp, rf, f"合計：同じ B を A+（サーバー側上限）と比べると低い {T['Bfull_vs_Aplus']['share_lower']:.0%}、中央値比 {T['Bfull_vs_Aplus']['median_ratio']:.1f}×", C_A, "集中 A+ の合計リスク /年"), ): ax.scatter(xa, yb, s=4, alpha=0.3, color=col) ax.plot([lo, hi], [lo, hi], "k--", lw=0.8, label="等リスク線"); ax.plot([lo, hi], [lo / 10, hi / 10], "k:", lw=0.8, label="分散側が 1/10") ax.set_xscale("log"); ax.set_yscale("log"); ax.set_xlim(1e-6, 1); ax.set_ylim(1e-7, 1) ax.set_xlabel(xl); ax.set_ylabel("分散側の合計リスク /年"); ax.set_title(title, fontsize=8); ax.legend(fontsize=7, loc="upper left") save(fig, "fig08_monte_carlo.png", f"図8　仮定表の LOW〜HIGH を三角分布で同時にサンプリングしたモンテカルロ（{summ['n_draws']:,} 回）。分散側の設計変数（100 ドメイン・権限上限 500 MW・e=0.3・ランプ 20%/分）は固定し、不確かな変数だけを振る。認証による p の低減は全構成に同じ倍率で適用する。左：制御面経路だけを見れば分散自律 B は常に低く、中央値で約 14 分の 1。左から2番目：ベンダー経路を合算すると、ソフトのみの B は A より低いものの中央値比 1.8 倍にとどまる。3番目：独立監視による下限保持とファームウェア系統の到達上限 12% を加えた B は A の約 9 分の 1。右：同じ B を、サーバー側に集計上限を置いた現実的な集中 A+ と比べると、低いのは約 8 割、中央値比は約 3 倍まで縮む。") # ---------------------------------------------------------------- fig09 tornado def fig09_tornado(): fig, axes = plt.subplots(1, 4, figsize=(19, 4.8)) jp = {"severity_k": "重大度の指数 k", "envelope_frac": "エンベロープ e", "ramp_limit_frac_per_min": "ランプ制限 /分", "p_compromise": "侵害確率 p", "remote_reach": "遠隔到達率", "grid_demand_mw": "系統需要", "frr_frac": "二次調整力 FRR", "fcr_frac": "一次調整力 FCR", "severity_x0": "重大度の尺度 x₀", "local_enforce": "ローカル強制生存率 ρ", "n_domains": "ドメイン数", "authority_cap_mw": "権限上限 A_max", "n_devices": "台数", "device_kw": "1台容量", "p_reduction_cert": "認証による p 低減", "t_recovery_h": "復旧時間", "sync_factor": "同期係数", "n_vendors": "ベンダー数", "vendor_top_share": "首位ベンダーシェア", "p_vendor_rel": "ベンダー経路の p 比", "vendor_hw_enforce": "独立監視の割合", "p_backend_rel": "バックエンド侵害の p 比", "domain_beta": "ドメイン間共通原因 β", "server_cap_mw": "サーバー側上限"} for ax, arch, col in ((axes[0], "A", C_A), (axes[1], "Aplus", C_A), (axes[2], "B", C_B), (axes[3], "Bfull", C_F)): t = rows(f"tornado_{arch}.csv")[:10] names = [jp.get(r["variable"], r["variable"]) for r in t][::-1] lo = np.array([fnum(r["low"]) for r in t])[::-1]; hi = np.array([fnum(r["high"]) for r in t])[::-1]; base = fnum(t[0]["base"]) y = np.arange(len(names)) ax.barh(y, hi - base, left=base, color=col, alpha=0.8, label="HIGH 値"); ax.barh(y, lo - base, left=base, color=col, alpha=0.35, label="LOW 値") ax.set_yticks(y); ax.set_yticklabels(names, fontsize=8); ax.axvline(base, color="k", lw=0.8) ax.set_xlabel("系統リスク /年"); ax.set_title(f"{ {'A': '集中 A', 'Aplus': '集中 A+（サーバー側上限）', 'B': '分散自律 B（ソフトのみ）', 'Bfull': 'B＋独立監視＋FW到達12%'}[arch] }：一変数感度（基準 {base:.2e}）", fontsize=8.5); ax.legend(fontsize=7) save(fig, "fig09_tornado.png", "図9　一変数感度（トルネード、フリート全体の合計リスク、設計変数は固定）。A は侵害確率と認証の低減倍率でほぼ決まり、容量側を動かしても重大度が飽和していて変わらない。A+ ではサーバー側上限とバックエンド侵害の確率比が加わる。ソフトのみの B は侵害確率・ベンダー経路の確率比・台数の順に効き、ドメイン分割は合計にはほとんど効かない。独立監視と到達上限を加えた B では、エンベロープ・系統規模・重大度の形（k, x₀）が効くようになる。どの構成でも侵害確率 p の幅が最大の不確かさである。") # ---------------------------------------------------------------- fig10 envelope sweep + vendor share limit def fig10_envelope_vendor(): es = rows("envelope_sweep.csv"); vs = rows("vendor_share_limit.csv") fig, axes = plt.subplots(1, 2, figsize=(13, 4.6)) ax = axes[0] for rho, col in ((0.8, "#bbb"), (0.9, C_H), (0.95, C_B), (0.99, "k")): pts = sorted((fnum(r["envelope_frac"]), fnum(r["x"])) for r in es if abs(fnum(r["local_enforce"]) - rho) < 1e-9) if pts: ax.plot([p[0] for p in pts], [p[1] for p in pts], marker="o", ms=3, color=col, label=f"ρ={rho}") ax.axhline(1, color="k", ls="--", lw=0.8); ax.text(0.02, 1.08, "x=1（FCR）", fontsize=7); ax.axhline(2.2, color=C_A, ls=":", lw=0.8); ax.text(0.02, 2.3, "x≈2.2（UFLS 開始）", fontsize=7, color=C_A) ax.set_xlabel("エンベロープ e（上位命令が動かせる割合）"); ax.set_ylabel("x（単一制御面を侵害されたとき）"); ax.set_yscale("log") ax.set_title("単一制御面のままエンベロープだけ付けた場合", fontsize=9); ax.legend(fontsize=8) ax = axes[1] nm = {"software": "ソフトのみ\n（下限もランプも FW と消える）", "independent": "独立監視で下限保持\n（ランプは消える）", "hypothetical": "仮想：下限もランプも残る\n（HW ランプ制限は未実在）"} names = [next(v for k, v in nm.items() if r["case"].startswith(k)) for r in vs] v1 = [fnum(r["max_vendor_share_for_x_below_1_0"]) * 100 for r in vs]; v15 = [fnum(r["max_vendor_share_for_x_below_1_5"]) * 100 for r in vs] xx = np.arange(len(names)); w = 0.35 ax.bar(xx - w / 2, v1, w, color=C_B, label="x<1（調整力内）"); ax.bar(xx + w / 2, v15, w, color=C_H, label="x<1.5（UFLS 前）") for i in range(len(names)): ax.annotate(f"{v1[i]:.1f}%", (xx[i] - w / 2, v1[i]), ha="center", va="bottom", fontsize=8); ax.annotate(f"{v15[i]:.1f}%", (xx[i] + w / 2, v15[i]), ha="center", va="bottom", fontsize=8) ax.set_xticks(xx); ax.set_xticklabels(names, fontsize=8); ax.set_ylabel("1 ファームウェア系統の到達上限（フリート比 %）"); ax.set_ylim(0, 50) ax.set_title("ファームウェア系統（鍵・OTA）の到達上限（基準：50 GW フリート・FCR 1,350 MW・FRR 2,250 MW）", fontsize=9); ax.legend(fontsize=8) save(fig, "fig10_envelope_vendor_limit.png", "図10　左：制御面を分割せずにエンベロープ（上位命令が動かせる割合 e）とローカル強制の生存率 ρ だけで x を下げようとした場合。ρ=0.95・e=0.3 でも x≈3.3 で UFLS 域に残り、エンベロープ単独では不十分。右：1 つのファームウェア系統（同じ署名鍵・同じ OTA 経路）の侵害が調整力を超えないための到達上限（§6 の x 定義：ステップ部分は FCR、全体は FCR+FRR と比較）。ソフトのみなら約 4%、独立監視で出力下限が残る設計なら約 11.5%、下限もランプも残る仮想設計なら約 31%。") # ---------------------------------------------------------------- fig11 graceful degradation ladder def fig11_ladder(): fig, ax = plt.subplots(figsize=(11, 5.2)); ax.set_xlim(0, 10); ax.set_ylim(0, 10); ax.axis("off"); ax.grid(False) steps = [ ("L0 正常", "上位の最適化命令を受理。市場・需給・混雑の最適化層が働く。", "#d9e8ee"), ("L1 上位を疑う", "命令が局所観測（V/f/RoCoF）と矛盾、またはエンベロープ超過 → 修正して受理。上位へ逸脱を報告。", "#cfe0e8"), ("L2 上位を無視", "命令が連続して拒否域、または認証失敗・証明書失効 → 直近の有効スケジュールと局所制御（droop・GFM）のみで運転。", "#c2d6e0"), ("L3 通信断", "上位との通信なし → 自律運転。既定の出力制御スケジュールまたはフェイルセーフ曲線。アイランド運転が可能ならセル単位で維持。", "#b3cad6"), ("L4 局所保護", "周波数・電圧が保護域 → ソフトに関係なくリレー／ハードウェア制限が動作。", "#a3bdcb"), ("L5 復旧", "正常化後、段階的に上位命令を再受理（ランプ制限付き）。再接続の同時性を抑える。", "#d9e8ee"), ] for i, (t, d, c) in enumerate(steps): y = 8.6 - i * 1.5 box(ax, 0.3, y, 1.9, 1.2, t, c, fs=9) ax.text(2.5, y + 0.6, d, va="center", fontsize=8.5) if i < len(steps) - 1: arrow(ax, 1.25, y, 1.25, y - 0.3) ax.text(5, 0.1, "各段で「残る機能」を定義しておくことが Graceful Degradation。L4 だけはソフトウェアの権限に依存しない。", ha="center", fontsize=8.5, color=C_B) save(fig, "fig11_graceful_degradation.png", "図11　劣化の階段（Graceful Degradation ladder）。上位の命令をどこまで信じるかを段階化し、各段で DER に残る機能を定義する。最終段 L4 はソフトウェアの権限に依存しないハードウェア保護で、ここが残る限り「侵害＝広域障害」にはならない。") # ---------------------------------------------------------------- fig12 incidents potential vs realised def fig12_incidents(): inc = rows("../data/incidents.csv") if (ROOT / "data" / "incidents.csv").exists() else [] if not inc: return inc = [r for r in inc if fnum(r["potential_mw"]) > 0] fig, ax = plt.subplots(figsize=(9.5, 5)) for r in inc: p, a = fnum(r["potential_mw"]), max(fnum(r["realised_mw"]), 0.5) col = C_A if r["kind"] == "cyber" else C_G ax.scatter(p, a, s=55, color=col, zorder=5, marker="o" if r["kind"] == "cyber" else "s") off = {"I01": (-8, 10), "I02": (-8, -14), "R01": (8, -12), "R02": (8, 6), "R05": (-110, -14), "R04": (8, 0), "R03": (8, 0), "I05": (8, -10)}.get(r["id"], (5, 4)) ax.annotate(f"{r['id']} {r['short_ja']}", (p, a), textcoords="offset points", xytext=off, fontsize=7, ha="right" if off[0] < 0 else "left") xx = np.logspace(-1, 6, 10); ax.plot(xx, xx, "k--", lw=0.8); ax.text(2e5, 1.3e5, "潜在＝顕在", fontsize=7, rotation=38) ax.axhline(1350, color=C_B, lw=0.8, ls=":"); ax.text(0.15, 1500, "基準系統の FCR", fontsize=7, color=C_B) ax.set_xscale("log"); ax.set_yscale("log"); ax.set_xlabel("潜在的に動かせた容量 MW（到達可能な台数×容量）"); ax.set_ylabel("実際に失われた・操作された容量 MW") ax.set_title("事例の潜在 Blast Radius と顕在化した影響：● サイバー　■ 非サイバーの参照事故", fontsize=9) save(fig, "fig12_incidents_potential_vs_realised.png", "図12　公開事例の「潜在的に動かせた容量」と「実際に動いた容量」。サイバー事例（●）は潜在 GW 級に対して実害がほぼゼロで、横軸と縦軸の乖離が今日の状況を表す。これは「安全だった」証拠ではなく、脆弱性が公表前に修正されたこと、そして分散電源のクラウドを経由した系統攻撃がまだ観測されていないことを示す（電力系統そのものへの攻撃は既にある）。非サイバーの参照事故（■）は同規模の擾乱が実際に起きたときの帰結を示す。") CAPTIONS_EN = { "fig01_architectures.png": "Figure 1. The three architectures compared. A: all DER on one national control plane; commands execute without device-side validation. H: ten regional domains, no device constraints. B: 100 independent credential domains, each device holding an output floor, ramp limit and authority cap, with the upper link advisory. Firmware distribution (OTA) is shared per vendor in all three.", "fig02_cbr_decomposition.png": "Figure 2. Decomposition of Cyber Blast Radius by compromise unit for the reference fleet (10 million × 5 kW). Installed (raw) capacity versus effective capacity after reach, authority cap and device floor. The national plane reaches 35 GW; the largest firmware line 10.5 GW regardless of domain splitting.", "fig03_severity_calibration.png": "Figure 3. Severity function S(x) = 1 − exp(−(x/1.7)³) fitted to literature points (black) with the swing-equation classification of the reference grid as blue steps. S is the probability of under-frequency load shedding or worse. Below x ≈ 0.5 the function is unidentified; the grey band shows k = 1–4.", "fig04_architecture_comparison.png": "Figure 4. Fourteen configurations, fleet-wide per year. Left: largest capacity moved by one compromise (solid = control plane / one domain; hatched = largest vendor firmware line). Middle: annual probability of a wide-area event (UFLS or worse) split into the control-plane channel (counting n domains and common cause β) and the vendor channel (cloud/OTA/firmware). Right: total. Domain splitting and the floor cut the plane channel by one to two orders, but the vendor channel bypasses authority caps and a software-only floor dies with the firmware, so software-only B totals about the same as certified A. What lowers the total is a floor held by an independent monitor (B5) and a firmware-line reach cap (B6, B7); a server-side aggregate cap (A+) brings the centralized design to the same order.", "fig05_attack_traces.png": "Figure 5. The same attack — a compromise of the control cloud or a firmware line cuts reachable DER output — applied to a single-area swing-equation model. Left: first 60 s of representative scenarios. Right: affected capacity versus frequency nadir for all scenarios. A1, A3, H1, B2, B5a and B5b reach load shedding; F1, B1 and B4 stay within reserves. In B5b the 30% floor survives but the ramp limit dies with the firmware, so 3.5 GW drops as a step. B5c (firmware-line reach capped at 12%) and A4 (server-side cap 1,500 MW) both land in the 'degraded' band just above the reserve.", "fig06_break_even.png": "Figure 6. Break-even between probability reduction and Blast Radius multiplication for a single compromise. Left: risk ratio A/B over the grid of probability reduction and Blast Radius multiplier. Right: the probability reduction a centralized design needs to equal a 10 MW distributed Blast Radius, for k = 1–4. Even at k = 1, 10,000 MW needs a 227-fold reduction.", "fig07_correlated_failure.png": "Figure 7. Monte Carlo of common-cause failure (vendor cloud, firmware line, aggregator credentials shared by groups of devices). Left: annual exceedance of simultaneous affected capacity. A single vendor (red) is rarely hit but its tail reaches 35 GW. Right: more vendors raise the probability that some vendor is compromised while lowering the 99th-percentile capacity. Diversity limits the size of an event, not its frequency; only with per-line reach caps and a surviving floor does it stay within reserves.", "fig08_monte_carlo.png": "Figure 8. Monte Carlo over the LOW–HIGH ranges of the uncertain assumptions (20,000 draws), with B's design variables fixed (100 domains, 500 MW cap, e = 0.3, 20%/min) and the certification factor applied to all configurations. Left: on the control-plane channel alone, distributed B is always lower (median 1/14). Second: adding the vendor channel, software-only B is lower than A but only by a median 1.8×. Third: B with an independent-monitor floor and a 12% firmware-line reach cap is 1/9 of A. Right: compared with a realistic centralized A+ (server-side aggregate cap), the same B is lower in about 83% of draws, median ratio 3.3.", "fig09_tornado.png": "Figure 9. One-variable sensitivity (tornado) of total fleet risk with design variables fixed. A is set almost entirely by compromise probability and the certification factor; capacity variables do not move it because severity is saturated. A+ adds the server-side cap and the backend-compromise ratio. Software-only B responds to p, the vendor-channel ratio and unit count; domain splitting barely affects the total. With the independent monitor and reach cap, the floor, grid size and severity shape (k, x₀) become influential. The range of p is the largest uncertainty everywhere.", "fig10_envelope_vendor_limit.png": "Figure 10. Left: trying to lower x with only a device floor e and local-enforcement survival ρ while keeping a single control plane; even at ρ = 0.95, e = 0.3, x ≈ 3.3 remains in the UFLS band. Right: the reach cap for one firmware line (same signing key, same OTA path) so that its compromise stays within reserves, using the Section 6 definition of x (step part against FCR, total against FCR+FRR): about 4% software-only, about 11.5% if an independent monitor keeps the floor, about 31% in a hypothetical design where floor and ramp both survive.", "fig11_graceful_degradation.png": "Figure 11. The graceful-degradation ladder: how far a device trusts upper-level commands at each level and what remains. L0 normal; L1 doubt (modify and report); L2 ignore (last valid schedule and local control); L3 link lost (autonomous, islanding where possible); L4 local protection (relay trips; independent monitor overrides below the floor); L5 staged recovery with ramp limits.", "fig12_incidents_potential_vs_realised.png": "Figure 12. Published incidents: capacity that could potentially have been moved versus capacity actually moved. Cyber incidents (●) show GW-scale potential with near-zero realised harm — evidence that vulnerabilities were fixed before publication and that no grid attack through a DER cloud has yet been observed, not evidence of safety; attacks on grids themselves already exist. Non-cyber reference events (■) show the consequence when a disturbance of that size actually occurs.", } def main() -> None: fig01_architectures(); fig02_cbr_decomposition(); fig03_severity(); fig04_architecture_comparison(); fig05_attack_traces() fig06_break_even(); fig07_ccf(); fig08_monte_carlo(); fig09_tornado(); fig10_envelope_vendor(); fig11_ladder(); fig12_incidents() (FIG / "captions.md").write_text("# 図版キャプション\n\n" + "\n\n".join(f"## {n}\n\n{c}" for n, c in captions) + "\n", encoding="utf-8") (FIG / "captions.en.md").write_text("# Figure captions\n\n" + "\n\n".join(f"## {n}\n\n{CAPTIONS_EN.get(n, c)}" for n, c in captions) + "\n", encoding="utf-8") print("captions written") if __name__ == "__main__": main()