#!/usr/bin/env python3 """ make_figures.py — 記事用の図を results/figures/ に生成する。 python3 model/make_figures.py すべて dwpt_model.py / monte_carlo.py の出力（results/*.csv, *.json）だけから描く。 図の数値を変えたいときは data/assumptions.csv を直して両スクリプトを再実行する。 """ from __future__ import annotations import csv import json import pathlib import sys import numpy as np import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt # noqa: E402 from matplotlib.ticker import FuncFormatter # noqa: E402 sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent)) import dwpt_model as M # noqa: E402 FIG = M.RESULTS / "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_DWPT, C_ALT, C_MCS, C_BAT, C_GRAY = "#c0392b", "#2c7fb8", "#2c7fb8", "#41ab5d", "#7f8c8d" def read_csv(p: pathlib.Path): with open(p, newline="", encoding="utf-8") as f: return list(csv.DictReader(f)) def fnum(x): return float(x) if x not in ("", "inf", None) else float("inf") def save(fig, name, caption): fig.tight_layout() fig.savefig(FIG / name, bbox_inches="tight") plt.close(fig) with open(FIG / "captions.md", "a", encoding="utf-8") as f: f.write(f"- `{name}`: {caption}\n") print("saved", name) def yen_fmt(x, _): if x >= 1e8: return f"{x/1e8:.0f}億" if x >= 1e4: return f"{x/1e4:.0f}万" return f"{x:.0f}" (FIG / "captions.md").write_text("# 図のキャプション（自動生成）\n\n", encoding="utf-8") base, low, high = M.load_assumptions("base"), M.load_assumptions("low"), M.load_assumptions("high") # ----------------------------------------------------------------------------- # 図1: 150 kW を kWh に変換する（通過時間 × 出力） # ----------------------------------------------------------------------------- rows = read_csv(M.RESULTS / "kw_to_kwh.csv") fig, ax = plt.subplots(figsize=(8, 4.6)) x = np.arange(len(rows)) kwh = [fnum(r["kwh_per_pass_duty100"]) for r in rows] kwh75 = [fnum(r["kwh_per_pass_duty75"]) for r in rows] ax.bar(x - 0.2, kwh, 0.4, color=C_DWPT, label="150 kW を区間全長で一定に受電（上限）") ax.bar(x + 0.2, kwh75, 0.4, color="#e59866", label="実効 75 %（位置ずれ・切替・SOC 制約）") for i, r in enumerate(rows): ax.annotate(f"{fnum(r['pass_time_s']):.0f} 秒\n{kwh[i]:.2f} kWh", (i - 0.2, kwh[i]), ha="center", va="bottom", fontsize=8) ax.set_xticks(x) ax.set_xticklabels([f"{int(fnum(r['length_m']))//1000 if fnum(r['length_m'])>=1000 else fnum(r['length_m'])/1000:g} km" if fnum(r["length_m"]) >= 1000 else f"{fnum(r['length_m']):.0f} m" for r in rows]) ax.set_yscale("log") ax.set_ylabel("1 回の通過で受け取る電力量 [kWh]（対数軸）") ax.set_xlabel("給電区間の長さ（100 km/h で通過）") ax.set_title("図1 「150 kW」を kWh に直す：300 m なら 10.8 秒・0.45 kWh（大型車の走行 0.3 km 分）") ax.legend(loc="upper left", fontsize=8) save(fig, "fig01_kw_to_kwh.png", "出力 150 kW は通過時間を掛けて初めて電力量になる。館山道 300 m 区間では 1 台あたり最大 0.45 kWh。") # ----------------------------------------------------------------------------- # 図2: 対象車への絞り込み（ファネル） # ----------------------------------------------------------------------------- fig, ax = plt.subplots(figsize=(9, 4.8)) labels = ["全交通量\n(AADT 上下計)", "大型車", "うち BEV", "うち DWPT\n受電器搭載", "うち敷設車線\nを走る車"] for j, (y, pen, col) in enumerate([(2026, 1.0, C_GRAY), (2030, 1.0, "#f39c12"), (2040, 0.50, "#8e44ad"), (2040, 1.0, C_DWPT)]): a = base vals = [a["aadt_total"], a["aadt_total"] * a["heavy_vehicle_fraction"]] vals.append(vals[-1] * a[f"ev_fraction_heavy_{y}"]) vals.append(vals[-1] * pen) vals.append(vals[-1] * a["lane_fraction"]) ax.plot(labels, vals, marker="o", color=col, label=f"{y} 年・受電器搭載率 {pen:.0%}（BASE）") ax.annotate(f"{vals[-1]:,.1f} 台/日", (4, vals[-1]), textcoords="offset points", xytext=(8, 0), fontsize=8, color=col) ax.set_yscale("log") ax.set_ylabel("台/日（対数軸）") ax.set_title(f"図2 館山道 君津PA付近：{base['aadt_total']:,.0f} 台/日（R3センサス平日）から DWPT 受電車は何台残るか") ax.legend(fontsize=8) save(fig, "fig02_target_funnel.png", "AADT から大型車・BEV・受電器搭載・敷設車線の順で絞ると、受電する車は 2040 年でも数百台/日の桁になる。") # ----------------------------------------------------------------------------- # 図3: コイル上在線時間（分/日）と年間供給量（kWh/m/年） # ----------------------------------------------------------------------------- grid = read_csv(M.RESULTS / "tateyama_grid_base.csv") fig, axes = plt.subplots(1, 2, figsize=(11, 4.6)) for y, col in [(2026, C_GRAY), (2030, "#f39c12"), (2035, "#8e44ad"), (2040, C_DWPT)]: rs = [r for r in grid if int(r["year"]) == y] pens = [fnum(r["dwpt_penetration"]) * 100 for r in rs] axes[0].plot(pens, [fnum(r["coil_minutes_per_day"]) for r in rs], marker="o", color=col, label=f"{y} 年") axes[1].plot(pens, [fnum(r["kwh_per_m_year"]) for r in rs], marker="o", color=col, label=f"{y} 年") axes[0].axhline(1440, color="k", ls="--", lw=0.8) axes[0].text(1.2, 1440 * 1.15, "1 日 = 1,440 分（コイルが常時使われる上限）", fontsize=8) for ax in axes: ax.set_xscale("log"); ax.set_yscale("log"); ax.set_xlabel("DWPT 受電器搭載率 [%]（BEV 大型車のうち）") ax.set_xticks([1, 5, 10, 25, 50, 100]); ax.get_xaxis().set_major_formatter(FuncFormatter(lambda v, _: f"{v:g}")) axes[0].set_ylabel("対象車がコイル上に存在する合計時間 [分/日]") axes[0].set_title("図3a 300 m のコイルに対象車がいる時間") axes[1].set_ylabel("年間供給電力量 [kWh/m/年]") axes[1].set_title("図3b 設備 1 m あたり年間供給量") axes[0].legend(fontsize=8) _g40 = next(r for r in grid if int(r["year"]) == 2040 and fnum(r["dwpt_penetration"]) == 1.0) save(fig, "fig03_coil_time_and_kwh_per_m.png", f"BASE 前提。2040 年・搭載率 100 % でも在線は 1 日 {fnum(_g40['coil_minutes_per_day']):.0f} 分（1,440 分の {fnum(_g40['coil_minutes_per_day'])/1440:.1%}）、年間 {fnum(_g40['kwh_per_m_year']):.0f} kWh/m。") # ----------------------------------------------------------------------------- # 図4: LCOS 対 対象車流量（損益分岐） # ----------------------------------------------------------------------------- flows = np.logspace(0, 4.5, 200) fig, ax = plt.subplots(figsize=(9, 5.2)) e1 = M.energy_per_pass_kwh(base["peak_power_kw"], 1000.0, base["speed_kmh"], base["coupling_duty"]) for name, cap, col, ls in [(f"LOW {low['capex_dwpt_yen_per_lane_km']/1e8:.1f} 億円/lane-km", low["capex_dwpt_yen_per_lane_km"], "#e59866", "--"), (f"BASE {base['capex_dwpt_yen_per_lane_km']/1e8:.1f} 億円", base["capex_dwpt_yen_per_lane_km"], C_DWPT, "-"), (f"HIGH {high['capex_dwpt_yen_per_lane_km']/1e8:.1f} 億円", high["capex_dwpt_yen_per_lane_km"], "#7b241c", ":")]: tot = [] for f in flows: l = M.dwpt_lcos_from(base, f, capex=cap) tot.append(l.total_yen_per_kwh) ax.plot(flows, tot, color=col, ls=ls, lw=2, label=f"DWPT 総 LCOS（CAPEX {name}）") alt26 = M.alternative_lcos(base, 2026)["total_yen_per_kwh"] alt40 = M.alternative_lcos(base, 2040)["total_yen_per_kwh"] ax.axhspan(alt40, alt26, color=C_ALT, alpha=0.15) ax.axhline(alt26, color=C_ALT, lw=1.5, label=f"代替：電池 +200 kWh ＋ MCS（2026: {alt26:.0f} 円/kWh）") ax.axhline(alt40, color=C_ALT, lw=1.5, ls="--", label=f"同 2040 年電池価格（{alt40:.0f} 円/kWh）") for y, pen, mk in [(2026, 1.0, "s"), (2030, 1.0, "^"), (2040, 0.5, "o"), (2040, 1.0, "D")]: f = M.target_flow_per_day(base["aadt_total"], base["heavy_vehicle_fraction"], base[f"ev_fraction_heavy_{y}"], pen, base["lane_fraction"]) l = M.dwpt_lcos_from(base, f) ax.plot([f], [l.total_yen_per_kwh], marker=mk, color="k", ms=8, ls="none") ax.annotate(f"館山道 {y} 年・搭載率 {pen:.0%}\n{f:,.0f} 台/日 → {l.total_yen_per_kwh:,.0f} 円/kWh", (f, l.total_yen_per_kwh), textcoords="offset points", xytext=(10, 6), fontsize=7.5) ax.set_xscale("log"); ax.set_yscale("log") ax.set_xlabel("DWPT で受電する車両数 [台/日・敷設車線]（対数軸）") ax.set_ylabel("車に届く 1 kWh あたりの費用 [円/kWh]（対数軸）") ax.set_title("図4 損益分岐：DWPT の LCOS が「電池＋MCS」に並ぶには何台/日が要るか") ax.legend(fontsize=7.5, loc="upper right") _be = {r["case"]: fnum(r["break_even_target_vehicles_per_day_per_lane"]) for r in read_csv(M.RESULTS / "break_even.csv") if int(r["year"]) == 2040} _f40 = M.target_flow_per_day(base["aadt_total"], base["heavy_vehicle_fraction"], base["ev_fraction_heavy_2040"], 1.0, base["lane_fraction"]) save(fig, "fig04_lcos_vs_flow.png", f"DWPT の LCOS は対象車流量に反比例する。2040 年の代替案 {alt40:.0f} 円/kWh に並ぶ流量は CAPEX LOW {_be['low_capex']:,.0f}／BASE {_be['base']:,.0f}／HIGH {_be['high_capex']:,.0f} 台/日・車線。館山道 2040 年・搭載率 100 % は {_f40:,.0f} 台/日で、BASE 損益分岐の {_f40/_be['base']:.0%}。") # ----------------------------------------------------------------------------- # 図5: 速度を落とすほど受電量が増える → 最適化すると「走行中」が消える # ----------------------------------------------------------------------------- rows = read_csv(M.RESULTS / "speed_sweep.csv") fig, ax = plt.subplots(figsize=(9, 4.4)) labels = [r["mode"].replace("stop", "停止") for r in rows] vals = [fnum(r["kwh"]) for r in rows] cols = [C_DWPT if int(float(r["speed_kmh"])) >= 50 else ("#f39c12" if int(float(r["speed_kmh"])) > 0 else C_ALT) for r in rows] ax.bar(labels, vals, color=cols) for i, v in enumerate(vals): ax.text(i, v * 1.15, f"{v:.2f}", ha="center", fontsize=8) ax.set_yscale("log") ax.set_ylabel("同じ 100 m・150 kW で受け取る kWh（対数軸）") ax.set_title("図5 受電量は速度に反比例する。設備利用率を上げる最適解は「止まること」") ax.text(0.01, 0.95, "赤: 高速本線 橙: ランプ・料金所・バス停 青: 停止（=静止充電）", transform=ax.transAxes, fontsize=8, va="top") save(fig, "fig05_speed_sweep.png", "同じ 100 m のコイルでも 100 km/h では 0.15 kWh、5 km/h で 3 kWh、30 分停止で 75 kWh。交通量の少ない区間で利用率を上げようとすると、コイルは停止点へ寄っていく。") # ----------------------------------------------------------------------------- # 図6: 区間を伸ばしても kWh/m/年と LCOS は変わらない（スケーリング則） # ----------------------------------------------------------------------------- rows = read_csv(M.RESULTS / "length_scaling.csv") fig, ax1 = plt.subplots(figsize=(8.5, 4.4)) L = [fnum(r["length_m"]) / 1000 for r in rows] ax1.plot(L, [fnum(r["kwh_per_pass"]) for r in rows], marker="o", color=C_DWPT, label="1 台あたり受電 kWh/pass（∝ 長さ）") ax1.plot(L, [fnum(r["capex_yen"]) / 1e8 for r in rows], marker="s", color=C_GRAY, label="CAPEX [億円]（∝ 長さ）") ax1.set_xscale("log"); ax1.set_yscale("log"); ax1.set_xlabel("給電区間長 [km]") ax1.set_ylabel("kWh/pass, 億円（対数軸）") ax2 = ax1.twinx() ax2.plot(L, [fnum(r["infra_yen_per_kwh"]) for r in rows], marker="D", color=C_ALT, lw=2.5, label="インフラ LCOS [円/kWh]（比例モデル：不変）") ax2.plot(L, [fnum(r["infra_yen_per_kwh_with_fixed_cost_and_cap"]) for r in rows], marker="v", color="#c0392b", lw=2, ls="--", label=f"同・区間固定費 {base['fixed_capex_per_segment_yen']/1e8:.1f} 億円＋受電上限 {base['usable_kwh_per_km']:.2f} kWh/km を加味") ax2.set_ylabel("円/kWh"); ax2.grid(False) ax2.set_ylim(0, max(fnum(r["infra_yen_per_kwh_with_fixed_cost_and_cap"]) for r in rows) * 1.3) h1, l1 = ax1.get_legend_handles_labels(); h2, l2 = ax2.get_legend_handles_labels() ax1.legend(h1 + h2, l1 + l2, fontsize=8, loc="upper left") ax1.set_title("図6 「道路を長くすれば解決」は成り立たない：受電量も費用も長さに比例し、円/kWh は動かない") _ls = {int(fnum(r["length_m"])): r for r in rows} save(fig, "fig06_length_scaling.png", f"2040 年・搭載率 100 %・BASE。費用も供給量も長さに比例するので、比例モデルでは区間長を 300 m から 30 km にしても円/kWh は {fnum(_ls[1000]['infra_yen_per_kwh']):.0f} 円のまま。受変電など区間固定費を入れると 300 m は {fnum(_ls[300]['infra_yen_per_kwh_with_fixed_cost_and_cap']):.0f} 円と割高になり、長くしても {fnum(_ls[30000]['infra_yen_per_kwh_with_fixed_cost_and_cap']):.0f} 円で下げ止まる。") # ----------------------------------------------------------------------------- # 図7: 電池を積む費用（円/kWh スループット） # ----------------------------------------------------------------------------- rows = read_csv(M.RESULTS / "battery_cases.csv") fig, ax = plt.subplots(figsize=(8.5, 4.4)) for y, col in [(2026, C_GRAY), (2030, "#f39c12"), (2035, "#8e44ad"), (2040, C_BAT)]: rs = [r for r in rows if int(r["year"]) == y] ax.plot([fnum(r["delta_kwh"]) for r in rs], [fnum(r["total_yen_per_kwh"]) for r in rs], marker="o", color=col, label=f"{y} 年電池価格") rs = [r for r in rows if int(r["year"]) == 2026] ax.fill_between([fnum(r["delta_kwh"]) for r in rs], 0, [fnum(r["throughput_yen_per_kwh"]) for r in rs], color=C_GRAY, alpha=0.15, label="うち電池償却分（2026）") ax.set_xlabel("追加する電池容量 ΔE [kWh]（基準 400 kWh に上乗せ）") ax.set_ylabel("電池 1 kWh を車に貯めて使う費用 [円/kWh]") _bv = [fnum(r["total_yen_per_kwh"]) for r in rows] ax.set_title(f"図7 「電池を積む」の費用：償却＋質量ペナルティ＋電費悪化でも {min(_bv):.0f}〜{max(_bv):.0f} 円/kWh") ax.legend(fontsize=8) save(fig, "fig07_battery_cases.png", "電池を +50〜+800 kWh 積んだときの 1 kWh スループット費用。DWPT が置き換えようとしている相手の価格。") # ----------------------------------------------------------------------------- # 図8: MCS LCOS 内訳（稼働率×重なり率） # ----------------------------------------------------------------------------- rows = read_csv(M.RESULTS / "mcs_lcos.csv") fig, ax = plt.subplots(figsize=(9, 4.6)) sel = [r for r in rows if fnum(r["overlap"]) in (0.0, 0.5, 1.0)] labels = [f"稼働率 {fnum(r['utilisation']):.0%}\n重なり {fnum(r['overlap']):.0%}" for r in sel] x = np.arange(len(sel)); bottom = np.zeros(len(sel)) for key, col, lab in [("infra_yen_per_kwh", C_MCS, "充電器＋系統連系 償却"), ("demand_charge_yen_per_kwh", "#74add1", "基本料金（kW）"), ("electricity_yen_per_kwh", "#abd9e9", "電力量料金"), ("time_yen_per_kwh", "#fdae61", "運転手・車両の待ち時間")]: v = np.array([fnum(r[key]) for r in sel]); ax.bar(x, v, bottom=bottom, color=col, label=lab); bottom += v for i, b in enumerate(bottom): ax.text(i, b + 2, f"{b:.0f}", ha="center", fontsize=8) ax.set_xticks(x); ax.set_xticklabels(labels, fontsize=7.5) ax.set_ylabel("円/kWh") ax.set_title("図8 MCS（1.2 MW）で 1 kWh 入れる費用の内訳：待ち時間は小さく、効くのは稼働率") ax.legend(fontsize=8) _mv = [fnum(r["total_yen_per_kwh"]) for r in rows]; _tv = [fnum(r["time_yen_per_kwh"]) for r in rows] save(fig, "fig08_mcs_lcos.png", f"MCS の LCOS は稼働率で {min(_mv):.0f}〜{max(_mv):.0f} 円/kWh。1 MW 級では時間費用は最大でも {max(_tv):.0f} 円/kWh で、重なり率はほとんど効かない。") # ----------------------------------------------------------------------------- # 図9: Dream Case からの歩み戻し（ウォーターフォール） # ----------------------------------------------------------------------------- d = M.dream_case_assumptions(low, base, high) def lcos_gap(a, year=2040): f = M.target_flow_per_day(a["aadt_total"], a["heavy_vehicle_fraction"], a[f"ev_fraction_heavy_{year}"], a["dwpt_penetration"], a["lane_fraction"]) l = M.dwpt_lcos_from(a, f) alt = M.alternative_lcos(a, year)["total_yen_per_kwh"] return l.total_yen_per_kwh, alt, f steps = [("Dream Case（全変数 DWPT 有利）", {})] cum = dict(d) for lab, upd in [ ("MCS 稼働率を BASE 15 % に", {"mcs_utilization": base["mcs_utilization"]}), ("MCS 費用を BASE に", {"mcs_charger_capex_yen": base["mcs_charger_capex_yen"], "mcs_site_grid_capex_yen": base["mcs_site_grid_capex_yen"]}), ("休憩との重なり 50 %", {"logistics_overlap": 0.5}), ("電池価格を 2040 BASE に", {"battery_pack_yen_per_kwh_2040": base["battery_pack_yen_per_kwh_2040"], "battery_cycle_life": base["battery_cycle_life"], "battery_kg_per_kwh": base["battery_kg_per_kwh"]}), ("敷設車線比率を BASE 45 %（片側 1 車線）", {"lane_fraction": base["lane_fraction"]}), ("受電器搭載率 50 %", {"dwpt_penetration": 0.5}), ("交通量・大型車率を実測 BASE に", {"aadt_total": base["aadt_total"], "heavy_vehicle_fraction": base["heavy_vehicle_fraction"]}), (f"BEV 比率 2040 BASE {base['ev_fraction_heavy_2040']:.0%}", {"ev_fraction_heavy_2040": base["ev_fraction_heavy_2040"]}), (f"速度 {base['speed_kmh']:.0f} km/h（大型貨物の規制速度）", {"speed_kmh": base["speed_kmh"]}), (f"CAPEX を BASE {base['capex_dwpt_yen_per_lane_km']/1e8:.1f} 億円", {"capex_dwpt_yen_per_lane_km": base["capex_dwpt_yen_per_lane_km"]}), ("寿命 20 年・O&M 2 %・割引率 5 %", {"lifetime_years": base["lifetime_years"], "om_fraction": base["om_fraction"], "discount_rate": base["discount_rate"]}), ]: steps.append((lab, upd)) labels, dw, al, fl = [], [], [], [] for lab, upd in steps: cum.update(upd) a_, b_, f_ = lcos_gap(cum) labels.append(lab); dw.append(a_); al.append(b_); fl.append(f_) fig, ax = plt.subplots(figsize=(10, 6)) yy = np.arange(len(labels)) ax.barh(yy - 0.18, dw, 0.36, color=C_DWPT, label="DWPT 総 LCOS") ax.barh(yy + 0.18, al, 0.36, color=C_ALT, label="電池＋MCS 総 LCOS") for i in range(len(labels)): ax.text(max(dw[i], al[i]) * 1.08, yy[i], f"DWPT {dw[i]:,.0f} vs 代替 {al[i]:,.0f} 円/kWh（対象 {fl[i]:,.0f} 台/日）", va="center", fontsize=7.5) ax.set_yticks(yy); ax.set_yticklabels(labels, fontsize=8); ax.invert_yaxis() ax.set_xscale("log"); ax.set_xlim(5, max(dw) * 40) ax.set_xlabel("円/kWh（対数軸）") ax.set_title("図9 Dream Case から変数を一つずつ現実へ戻す：どこで DWPT が負けに転じるか（2040 年）") ax.legend(fontsize=8, loc="lower right") _flip = next((labels[i] for i in range(len(labels)) if dw[i] > al[i]), None) save(fig, "fig09_dream_walkback.png", f"DWPT に極端に有利な前提では {dw[0]:.0f} 対 {al[0]:.0f} 円/kWh で勝つ。この順で変数を現実値へ戻すと「{_flip}」の段で逆転する。順序を変えると逆転する段は変わる（図9b の一変数感度を併読）。") with open(M.RESULTS / "dream_walkback.csv", "w", newline="", encoding="utf-8") as f: w = csv.writer(f); w.writerow(["step", "dwpt_total_yen_per_kwh", "alt_total_yen_per_kwh", "target_flow_per_day"]) for i in range(len(labels)): w.writerow([labels[i], f"{dw[i]:.1f}", f"{al[i]:.1f}", f"{fl[i]:.1f}"]) # ----------------------------------------------------------------------------- # 図9b: 一変数感度（トルネード）：Dream Case から 1 変数だけ BASE に戻したときの DWPT − 代替 の差 # ----------------------------------------------------------------------------- dw0, al0, f0 = lcos_gap(d) one_at_a_time = [] for lab, upd in steps[1:]: a_, b_, f_ = lcos_gap(dict(d, **upd)) one_at_a_time.append((lab, a_, b_, f_, a_ - b_)) one_at_a_time.sort(key=lambda t: -t[4]) fig, ax = plt.subplots(figsize=(10, 5.2)) yy = np.arange(len(one_at_a_time)) vals = [t[4] for t in one_at_a_time] ax.barh(yy, vals, color=[C_DWPT if v > 0 else C_ALT for v in vals]) ax.axvline(0, color="k", lw=1) for i, t in enumerate(one_at_a_time): ax.text(3, yy[i], f"DWPT {t[1]:,.0f} vs 代替 {t[2]:,.0f} 円/kWh（{t[3]:,.0f} 台/日）", va="center", ha="left", fontsize=7.5) ax.set_xlim(min(vals) * 1.05, abs(min(vals)) * 0.75) ax.set_yticks(yy); ax.set_yticklabels([t[0] for t in one_at_a_time], fontsize=8); ax.invert_yaxis() ax.set_xlabel(f"DWPT 総 LCOS − 代替 総 LCOS [円/kWh]（Dream Case 基準 {dw0 - al0:,.0f}。正なら DWPT が負け）") ax.set_title("図9b Dream Case から 1 変数だけ現実に戻す：単独で勝敗をひっくり返す変数はどれか（2040 年）") _single_flips = [t[0] for t in one_at_a_time if t[4] > 0] save(fig, "fig09b_tornado.png", f"他の全変数を DWPT 有利のままにしても、単独で逆転させる変数は {len(_single_flips)} 個：{'、'.join(_single_flips) if _single_flips else 'なし'}。") with open(M.RESULTS / "dream_tornado.csv", "w", newline="", encoding="utf-8") as f: w = csv.writer(f); w.writerow(["variable_reset_to_base", "dwpt_total_yen_per_kwh", "alt_total_yen_per_kwh", "target_flow_per_day", "dwpt_minus_alt"]) for t in one_at_a_time: w.writerow([t[0], f"{t[1]:.1f}", f"{t[2]:.1f}", f"{t[3]:.1f}", f"{t[4]:.1f}"]) # Dream の技術・金融前提のまま、交通量だけ館山道 BASE（レビュー R3-03） _d_tat = dict(d, aadt_total=base["aadt_total"], heavy_vehicle_fraction=base["heavy_vehicle_fraction"], ev_fraction_heavy_2040=base["ev_fraction_heavy_2040"], lane_fraction=base["lane_fraction"]) _dt = lcos_gap(_d_tat) _d_tat_altbase = dict(_d_tat, mcs_utilization=base["mcs_utilization"], mcs_charger_capex_yen=base["mcs_charger_capex_yen"], mcs_site_grid_capex_yen=base["mcs_site_grid_capex_yen"], battery_pack_yen_per_kwh_2040=base["battery_pack_yen_per_kwh_2040"], battery_cycle_life=base["battery_cycle_life"], battery_kg_per_kwh=base["battery_kg_per_kwh"], logistics_overlap=0.5) _dt2 = lcos_gap(_d_tat_altbase) with open(M.RESULTS / "dream_at_tateyama_flow.json", "w", encoding="utf-8") as f: json.dump({"dream_tech_finance_tateyama_base_traffic_pen100": {"dwpt_total_yen_per_kwh": round(_dt[0], 1), "alt_total_yen_per_kwh_dream_alt": round(_dt[1], 1), "target_flow_per_day": round(_dt[2], 1)}, "same_but_alternative_at_base": {"dwpt_total_yen_per_kwh": round(_dt2[0], 1), "alt_total_yen_per_kwh": round(_dt2[1], 1), "target_flow_per_day": round(_dt2[2], 1)}}, f, ensure_ascii=False, indent=2) # ----------------------------------------------------------------------------- # 図10: Monte Carlo — 勝率対流量、散布図 # ----------------------------------------------------------------------------- curve = read_csv(M.RESULTS / "monte_carlo_winrate_by_flow.csv") mc = read_csv(M.RESULTS / "monte_carlo.csv") summ = json.load(open(M.RESULTS / "monte_carlo_summary.json", encoding="utf-8")) fig, axes = plt.subplots(1, 2, figsize=(12, 4.8)) xs = [np.sqrt(fnum(c["flow_bin_low"]) * fnum(c["flow_bin_high"])) for c in curve] axes[0].plot(xs, [fnum(c["dwpt_win_rate"]) * 100 for c in curve], marker="o", color=C_DWPT) axes[0].set_xscale("log"); axes[0].set_xlabel("DWPT 受電車 [台/日・車線]"); axes[0].set_ylabel("DWPT が安くなる確率 [%]") axes[0].set_title(f"図10a Monte Carlo {summ['n_cases']:,} ケース：勝率は流量でほぼ決まる（全体勝率 {summ['dwpt_win_rate']*100:.2f} %）") f_ = np.array([fnum(r["target_flow_per_day"]) for r in mc]); c_ = np.array([fnum(r["capex_yen_per_lane_km"]) for r in mc]); w_ = np.array([r["winner"] == "DWPT" for r in mc]) axes[1].scatter(f_[~w_], c_[~w_] / 1e8, s=3, color=C_ALT, alpha=0.25, label="電池＋MCS が安い") axes[1].scatter(f_[w_], c_[w_] / 1e8, s=10, color=C_DWPT, alpha=0.9, label="DWPT が安い") axes[1].set_xscale("log"); axes[1].set_yscale("log") axes[1].set_xlabel("DWPT 受電車 [台/日・車線]"); axes[1].set_ylabel("DWPT CAPEX [億円/lane-km]") axes[1].set_title("図10b 勝ちケースの位置：高流量 × 低 CAPEX の隅にだけ存在する") axes[1].legend(fontsize=8, loc="lower left") save(fig, "fig10_monte_carlo.png", f"館山道型の交通分布で {summ['n_cases']:,} ケースを引くと DWPT が安くなる確率は {summ['dwpt_win_rate']*100:.1f} %。右図の赤点（勝ちケース）は高流量 × 低 CAPEX の隅にしか現れない。") # ----------------------------------------------------------------------------- # 図11: 年間費用と年間供給量の対比（館山道 300 m を 1 km に正規化） # ----------------------------------------------------------------------------- grid = read_csv(M.RESULTS / "tateyama_grid_base.csv") fig, ax = plt.subplots(figsize=(9, 4.6)) ann = M.annualised_cost(base["capex_dwpt_yen_per_lane_km"], base["discount_rate"], base["lifetime_years"], base["om_fraction"]) cases = [(2026, 1.0), (2030, 0.1), (2030, 1.0), (2040, 0.1), (2040, 0.5), (2040, 1.0)] labels = [f"{y}年\n搭載率{p:.0%}" for y, p in cases] rev = [] for y, p in cases: r = next(rr for rr in grid if int(rr["year"]) == y and abs(fnum(rr["dwpt_penetration"]) - p) < 1e-9) rev.append(fnum(r["mwh_per_km_year"]) * 1000 * M.alternative_lcos(base, y)["total_yen_per_kwh"]) x = np.arange(len(cases)) ax.bar(x - 0.2, [ann] * len(cases), 0.4, color=C_DWPT, label="DWPT 1 lane-km の年額費用（償却＋O&M, BASE）") ax.bar(x + 0.2, rev, 0.4, color=C_ALT, label="その 1 km が 1 年に供給する電力量を代替案単価で評価した価値") ax.set_xticks(x); ax.set_xticklabels(labels, fontsize=8) ax.set_yscale("log"); ax.yaxis.set_major_formatter(FuncFormatter(yen_fmt)) ax.set_ylabel("円/年（対数軸）") ax.set_title("図11 館山道で 1 km 敷設したら：年に払う額 vs 年に生む価値") ax.legend(fontsize=8) save(fig, "fig11_annual_cost_vs_value.png", f"BASE CAPEX の年額は約 {ann/1e4:,.0f} 万円/lane-km。2040 年・搭載率 100 % でも生む価値は約 {rev[-1]/1e4:,.0f} 万円/lane-km で届かない。2030 年以前は 2〜3 桁足りない。") # ----------------------------------------------------------------------------- # 図12: MCS ハブの電力集中（系統連系＋バッファ）と DWPT の分散電源の比較 # ----------------------------------------------------------------------------- rows = read_csv(M.RESULTS / "mcs_hub.csv") fig, ax = plt.subplots(figsize=(8.5, 4.4)) n = [int(float(r["simultaneous_trucks"])) for r in rows] ax.bar([str(v) for v in n], [fnum(r["peak_mw"]) for r in rows], color=C_MCS) for i, r in enumerate(rows): ax.text(i, fnum(r["peak_mw"]) + 1, f"{fnum(r['peak_mw']):.0f} MW\n連系＋1h電池 {fnum(r['total_yen'])/1e8:.0f} 億円\n→ {fnum(r['yen_per_kwh_delivered']):.1f} 円/kWh", ha="center", fontsize=8) ax.set_xlabel("同時充電台数（1.2 MW/台）") ax.set_ylabel("ピーク需要 [MW]") ax.set_ylim(0, max(fnum(r["peak_mw"]) for r in rows) * 1.5) ax.set_title("図12 SA/PA 電力集中問題：100 台同時でも 120 MW、連系＋バッファは供給 kWh あたり数円") save(fig, "fig12_mcs_hub.png", "MCS ハブの系統連系と電池バッファの費用は、供給電力量で割ると数円/kWh に収まる。DWPT 側の『分散できる』利点はこの数円が上限。") # ----------------------------------------------------------------------------- # 図13: 汎用回廊の DWPT 勝利空間（流量 × CAPEX ヒートマップ） # ----------------------------------------------------------------------------- wide_p = M.RESULTS / "monte_carlo_wide_grid.npz" if wide_p.exists(): z = np.load(wide_p) grid_, fb, cb = z["grid"], z["flow_bins"], z["capex_bins"] fig, ax = plt.subplots(figsize=(10, 5.6)) pc = ax.pcolormesh(fb, cb / 1e8, grid_ * 100, cmap="RdYlBu_r", vmin=0, vmax=100, shading="flat") cbar = fig.colorbar(pc, ax=ax); cbar.set_label("DWPT が「電池＋MCS」より安くなる確率 [%]") ax.set_xscale("log"); ax.set_yscale("log") ax.set_xlabel("DWPT 受電車 [台/日・敷設車線]（他の全変数は LOW〜HIGH で乱数）") ax.set_ylabel("DWPT CAPEX [億円/lane-km]") # 館山道の点と、新東名級専用レーンの点 tf = lambda y, pen: M.target_flow_per_day(base["aadt_total"], base["heavy_vehicle_fraction"], base[f"ev_fraction_heavy_{y}"], pen, base["lane_fraction"]) pts = [("館山道 2030・搭載率100%", tf(2030, 1.0), "s", (-10, 14)), ("館山道 2040・搭載率100%", tf(2040, 1.0), "D", (-10, -18)), ("館山道 2040・搭載率10%", tf(2040, 0.1), "o", (-80, -18)), ("新東名級 大型BEV 全車搭載・専用レーン（仮）", 17000, "^", (-150, 14))] for lab, f_, mk, off in pts: ax.plot([f_], [base["capex_dwpt_yen_per_lane_km"] / 1e8], marker=mk, color="k", ms=9, ls="none") ax.annotate(lab, (f_, base["capex_dwpt_yen_per_lane_km"] / 1e8), textcoords="offset points", xytext=off, fontsize=7.5) ax.axhline(low["capex_dwpt_yen_per_lane_km"] / 1e8, color="k", ls=":", lw=0.8); ax.text(1.2, low["capex_dwpt_yen_per_lane_km"] / 1e8 * 1.05, "文献最低 CAPEX", fontsize=7.5) ax.axhline(high["capex_dwpt_yen_per_lane_km"] / 1e8, color="k", ls=":", lw=0.8); ax.text(1.2, high["capex_dwpt_yen_per_lane_km"] / 1e8 * 1.05, "実証プロジェクト実績級 CAPEX", fontsize=7.5) ax.set_title("図13 DWPT が勝つパラメータ空間：右下（高流量 × 低 CAPEX）だけが赤い") save(fig, "fig13_win_space_heatmap.png", "汎用回廊 10 万ケース。DWPT が安くなる確率が 50 % を超えるのは、対象車 数千台/日・車線 かつ CAPEX 2 億円/lane-km 以下の領域に限られる。館山道の点はすべて青の領域にある。") # ----------------------------------------------------------------------------- # 図14: 海外実証の時系列（開始 → 終了/凍結/継続） # ----------------------------------------------------------------------------- # (名称, 開始年, 終了年 or None=継続, 状態, 方式) TIMELINE = [ ("韓国 KAIST OLEV（亀尾・世宗 ほか）", 2009, 2021, "撤去", "無線"), ("スウェーデン eRoadArlanda（導電レール）", 2018, 2021, "終了", "導電"), ("スウェーデン Smartroad Gotland（無線）", 2020, 2023, "終了・撤去", "無線"), ("スウェーデン E20 恒久ERS 計画（21 km）", 2021, 2025, "入札中止→国家計画から削除", "未定"), ("ドイツ eHighway A5 ELISA（架線）", 2019, 2024, "全停止", "架線"), ("ドイツ eHighway A1 FESH（架線）", 2019, 2024, "全停止", "架線"), ("ドイツ eWayBW B462（架線）", 2021, 2024, "停止→2025 撤去", "架線"), ("イタリア Arena del Futuro（A35 無線）", 2021, None, "デモ継続・公道化なし", "無線"), ("ノルウェー Trondheim バス（無線 82 m）", 2023, None, "次期契約で要件化せず", "無線"), ("米国 Detroit 14th St（無線 400 m）", 2023, 2025, "試験完了・延伸未着工", "無線"), ("米国 Indiana US-231（無線）", 2024, None, "最終報告・標準化へ", "無線"), ("フランス A10 Charge as you drive（無線）", 2025, None, "実証中", "無線"), ("日本 柏の葉 市道（無線 10 kW）", 2023, None, "継続", "無線"), ("日本 大阪・関西万博 EVバス（無線 30 kW）", 2025.3, 2025.8, "会期終了（184 日）", "無線"), ("日本 館山道 本線（無線 150 kW級・300 m）", 2027.3, 2029, "計画（2027 年度以降・複数回）", "無線"), ] fig, ax = plt.subplots(figsize=(10.5, 6.2)) colors = {"無線": C_DWPT, "架線": "#8e44ad", "導電": C_BAT, "未定": "#999999"} for i, (name, s, e, status, kind) in enumerate(reversed(TIMELINE)): end = e if e is not None else 2026.75 planned = s > 2026.75 ax.barh(i, end - s, left=s, color=colors.get(kind, "#999"), alpha=0.85 if (e and not planned) else 0.45, height=0.6, hatch="//" if planned else None) ax.text(end + 0.15, i, status, va="center", fontsize=7.5, color="#333") ax.set_yticks(range(len(TIMELINE))); ax.set_yticklabels([t[0] for t in reversed(TIMELINE)], fontsize=8) ax.axvline(2026.75, color="k", ls="--", lw=0.8); ax.text(2026.8, len(TIMELINE) - 0.4, "2026-10", fontsize=7.5) ax.set_xlim(2008, 2031); ax.set_xlabel("年") ax.set_title("図14 走行中給電（ERS）実証の時系列：終了・凍結・撤去が並び、商用公道区間は無い") from matplotlib.patches import Patch ax.legend(handles=[Patch(color=C_DWPT, label="無線"), Patch(color="#8e44ad", label="架線"), Patch(color=C_BAT, label="導電レール"), Patch(color="#999", label="方式未定")], loc="lower left", fontsize=8) save(fig, "fig14_international_timeline.png", "2009 年以降の主要 ERS 実証。濃い帯は終了済み、薄い帯は継続中。2026 年 10 月時点で、実証から商用運用に移った公道区間は確認できない。") # ----------------------------------------------------------------------------- # 図15: 円/vehicle-km と 円/tonne-km（館山道 2040 年、搭載率別） # ----------------------------------------------------------------------------- grid_rows = [r for r in read_csv(M.RESULTS / "tateyama_grid_base.csv") if r["year"] == "2040"] pens = [fnum(r["dwpt_penetration"]) for r in grid_rows] kwh_km = base["truck_kwh_per_km"] payload_t = 10.0 # 大型トラックの平均積載量（t）: 積載効率 ~40 % × 最大積載 ~13 t を丸めた値 dwpt_vkm = [fnum(r["lcos_total_yen_per_kwh"]) * kwh_km for r in grid_rows] alt = M.alternative_lcos(base, 2040)["total_yen_per_kwh"] * kwh_km fig, ax = plt.subplots(figsize=(8.5, 4.6)) x = np.arange(len(pens)) ax.bar(x - 0.2, dwpt_vkm, width=0.4, color=C_DWPT, label="DWPT（インフラ＋受電器＋電気）") ax.bar(x + 0.2, [alt] * len(pens), width=0.4, color=C_MCS, label="電池＋MCS") for i, v in enumerate(dwpt_vkm): ax.text(i - 0.2, v * 1.15, f"{v:,.0f} 円/km\n{v / payload_t:,.0f} 円/t-km", ha="center", fontsize=7.5) ax.text(len(pens) - 1 + 0.2, alt * 1.15, f"{alt:,.0f} 円/km\n{alt / payload_t:,.1f} 円/t-km", ha="center", fontsize=7.5) ax.axhline(25 * payload_t, color="k", ls=":", lw=0.8); ax.text(-0.4, 25 * payload_t * 1.1, f"トラック運賃 25 円/t-km × {payload_t:.0f} t = {25 * payload_t:,.0f} 円/km（エネルギー費はこの一部）", fontsize=7.5) ax.set_yscale("log"); ax.set_ylim(50, max(dwpt_vkm) * 4); ax.set_xticks(x); ax.set_xticklabels([f"搭載率 {p:.0%}" for p in pens]) ax.set_ylabel("エネルギー供給費 [円/vehicle-km]（対数）") ax.set_title(f"図15 館山道 2040 年：1 台が 1 km 走るためのエネルギー供給費（電費 {kwh_km} kWh/km）") ax.legend(fontsize=8, loc="upper right") save(fig, "fig15_yen_per_vehicle_km.png", f"円/kWh を電費 {kwh_km} kWh/km で車両 km に換算。搭載率 100 % でも DWPT は 1 km あたり約 {dwpt_vkm[-1]:,.0f} 円、電池＋MCS は約 {alt:,.0f} 円。運賃 25 円/t-km（10 t 積載で 250 円/km）と比べると、DWPT のエネルギー費は運賃そのものを上回る。") print("done →", FIG)