#!/usr/bin/env python3 """ dwpt_model.py — 高速道路 走行中無線給電（DWPT）の物理・経済モデル 目的 ---- 「150 kW」という宣伝数字を、 通過時間 [s] → 1台あたり受電量 [kWh/pass] → 年間供給量 [kWh/m/year] → 設備費の年額化 [円/km/year] → 供給電力量あたり単価 LCOS [円/kWh] へ変換し、同じ 1 kWh を車に届ける代替手段（電池を積む＋MCS 急速充電、架線式 ERS）と 同じ単位で比べる。すべての式は本ファイルに書き、`data/assumptions.csv` の LOW/BASE/HIGH を差し替えれば誰でも再計算できる。 単位の約束 ---------- 長さ m、速度 km/h（内部で m/s へ）、時間 s / h / year、電力 kW、電力量 kWh、金額 円。 交通量は「台/日（上下合計）」。年は 365 日、1 年 = 8,760 時間。 使い方 ------ python3 model/dwpt_model.py # 全ケースを計算し results/*.csv を出力 from dwpt_model import * # 関数単位で利用 """ from __future__ import annotations import csv import json import math import pathlib from dataclasses import dataclass, asdict, replace from typing import Dict, Iterable, List ROOT = pathlib.Path(__file__).resolve().parents[1] DATA = ROOT / "data" RESULTS = ROOT / "results" SEC_PER_HOUR = 3600.0 SEC_PER_DAY = 86400.0 DAYS_PER_YEAR = 365.0 HOURS_PER_YEAR = 8760.0 # 受電器搭載車が 1 年に走る DWPT 敷設区間の距離（回廊が整備された場合の DWPT 有利側の既定値） DWPT_KM_PER_VEHICLE_YEAR_CORRIDOR = 25000.0 # ============================================================================= # 0. 前提値の読み込み # ============================================================================= def load_assumptions(case: str = "base", path: pathlib.Path = DATA / "assumptions.csv") -> Dict[str, float]: """assumptions.csv から LOW/BASE/HIGH の 1 列を dict で返す。""" assert case in ("low", "base", "high") out: Dict[str, float] = {} with open(path, newline="", encoding="utf-8") as f: for row in csv.DictReader(f): try: out[row["parameter"]] = float(row[case]) except ValueError: pass return out # ============================================================================= # 1. 物理：時間・占有率・エネルギー # ============================================================================= def kmh_to_ms(v_kmh: float) -> float: return v_kmh / 3.6 def pass_time_s(section_m: float, v_kmh: float) -> float: """区間長 D [m] を速度 v で通過する時間 [s]。 t = D / v""" return section_m / kmh_to_ms(v_kmh) def occupancy_time_s(section_m: float, vehicle_len_m: float, v_kmh: float) -> float: """車体の一部がコイル区間のどこかに重なっている時間 [s]。 t = (D + L) / v （先頭がコイル始端に入ってから、末尾がコイル終端を抜けるまで）""" return (section_m + vehicle_len_m) / kmh_to_ms(v_kmh) def point_occupancy(flow_per_day: float, vehicle_len_m: float, v_kmh: float) -> float: """道路上の一点が車両に覆われている時間割合（無次元）。 occupancy ≈ q · L / v （q: 台/s, L: m, v: m/s） 例: 38,000 台/日・L=4.5 m・v=100 km/h → 0.44 台/s × 4.5 m / 27.8 m/s ≈ 0.071 (7.1 %)""" q_per_s = flow_per_day / SEC_PER_DAY return q_per_s * vehicle_len_m / kmh_to_ms(v_kmh) def energy_per_pass_kwh(power_kw: float, section_m: float, v_kmh: float, duty: float = 1.0) -> float: """1 回の通過で受け取る電力量 [kWh]。 E = P · duty · t / 3600 duty: ピーク出力に対する実効平均出力比（位置ずれ・セグメント切替・SOC 受入制限）。 例: 150 kW・300 m・100 km/h・duty 1.0 → 10.8 s → 0.45 kWh""" return power_kw * duty * pass_time_s(section_m, v_kmh) / SEC_PER_HOUR def target_flow_per_day(aadt: float, heavy_frac: float, ev_frac: float, dwpt_pen: float, lane_frac: float) -> float: """DWPT で実際に受電する車両数 [台/日]。 = AADT × 大型車比率 × BEV 比率 × DWPT 受電器搭載率 × 敷設車線を走る割合""" return aadt * heavy_frac * ev_frac * dwpt_pen * lane_frac def annual_energy_kwh(flow_per_day: float, e_pass_kwh: float) -> float: return flow_per_day * e_pass_kwh * DAYS_PER_YEAR def kwh_per_m_year(annual_kwh: float, section_m: float) -> float: """設備 1 m あたり年間供給量 [kWh/m/year]。本記事の主指標。""" return annual_kwh / section_m def mwh_per_km_year(annual_kwh: float, section_m: float) -> float: return kwh_per_m_year(annual_kwh, section_m) * 1000.0 / 1000.0 # kWh/m = MWh/km def coil_time_share(flow_per_day: float, section_m: float, vehicle_len_m: float, v_kmh: float) -> float: """対象車がコイル区間上に存在する時間の合計 ÷ 1 日。 1 日の合計秒 = flow × (D+L)/v。これを 86,400 s で割る。1 を超える場合は複数台同時在線。""" return flow_per_day * occupancy_time_s(section_m, vehicle_len_m, v_kmh) / SEC_PER_DAY def power_capacity_utilisation(annual_kwh: float, power_kw: float) -> float: """設備容量利用率 = 年間供給 kWh ÷ (ピーク kW × 8,760 h)。 300 m 区間では同時在線が稀なので、設備容量をピーク出力 1 台分とみなす。""" return annual_kwh / (power_kw * HOURS_PER_YEAR) # ============================================================================= # 2. 財務：年額化と LCOS # ============================================================================= def crf(rate: float, years: float) -> float: """資本回収係数 CRF = r / (1 − (1+r)^−n)。r=0 なら 1/n。""" if rate <= 0: return 1.0 / years return rate / (1.0 - (1.0 + rate) ** (-years)) def annualised_cost(capex: float, rate: float, years: float, om_fraction: float) -> float: """年額費用 = CAPEX × (CRF + O&M 比率)""" return capex * (crf(rate, years) + om_fraction) @dataclass class DwptLcos: infra_yen_per_kwh: float receiver_yen_per_kwh: float electricity_yen_per_kwh: float total_yen_per_kwh: float annual_kwh_per_lane_km: float annualised_capex_yen_per_lane_km: float demand_charge_yen_per_kwh: float = 0.0 contract_kw_per_lane_km: float = 0.0 def dwpt_contract_kw(flow_per_day: float, power_kw: float, v_kmh: float, min_kw: float, peak_hour_share: float = 0.10, length_m: float = 1000.0, duty: float = 0.9) -> float: """区間（1 lane-km）の契約電力 [kW]。 日本の高圧契約（実量制）は 30 分平均の最大需要電力（デマンド値）で契約電力が決まるので、 12 秒の 150 kW パルスはそのまま契約 kW にはならない。ピーク時間帯の 30 分間に通る受電車の 受電 kWh を 0.5 h で割った平均 kW を取り、下限 min_kw（高圧契約の最小 50 kW、または保守的に 1 台分）を適用する。 瞬時 150〜300 kW を受ける変圧器・配電設備の費用は fixed_capex_per_segment_yen（CAPEX 側）で数える。""" e_pass = energy_per_pass_kwh(power_kw, length_m, v_kmh, duty) peak_half_hour_vehicles = flow_per_day * peak_hour_share / 2.0 avg_kw = peak_half_hour_vehicles * e_pass / 0.5 return max(min_kw, avg_kw) def dwpt_lcos( capex_yen_per_lane_km: float, rate: float, lifetime_years: float, om_fraction: float, flow_per_day: float, power_kw: float, v_kmh: float, duty: float, eff: float, electricity_yen_per_kwh: float, receiver_cost_yen: float, receiver_life_years: float, receiver_kwh_per_year: float, demand_charge_yen_per_kw_month: float = 0.0, contract_kw_min: float = 150.0, peak_hour_share: float = 0.10, ) -> DwptLcos: """DWPT の供給電力量あたり単価 [円/kWh]。 - infra: 1 lane-km の年額費用 ÷ その 1 lane-km が 1 年に車へ渡す kWh （kWh/pass ∝ D、CAPEX ∝ D なので D に依存しない。1 km で計算） - receiver: 車載受電器の年額 ÷ その車が 1 年に DWPT で受け取る kWh - electricity: 購入電力 ÷ 効率 - demand charge: 契約電力（ピーク同時在線台数×kW、下限 1 台分）× 基本料金 ÷ 年間 kWh（実量制。demand_charge=0 なら未計上）""" e_pass_1km = energy_per_pass_kwh(power_kw, 1000.0, v_kmh, duty) annual_kwh = annual_energy_kwh(flow_per_day, e_pass_1km) ann_capex = annualised_cost(capex_yen_per_lane_km, rate, lifetime_years, om_fraction) infra = ann_capex / annual_kwh if annual_kwh > 0 else float("inf") recv = annualised_cost(receiver_cost_yen, rate, receiver_life_years, 0.0) / receiver_kwh_per_year if receiver_kwh_per_year > 0 else float("inf") elec = electricity_yen_per_kwh / eff contract_kw = dwpt_contract_kw(flow_per_day, power_kw, v_kmh, contract_kw_min, peak_hour_share, 1000.0, duty) if demand_charge_yen_per_kw_month > 0 else 0.0 demand = contract_kw * demand_charge_yen_per_kw_month * 12.0 / annual_kwh if (annual_kwh > 0 and demand_charge_yen_per_kw_month > 0) else 0.0 return DwptLcos(infra, recv, elec, infra + recv + elec + demand, annual_kwh, ann_capex, demand, contract_kw) def dwpt_lcos_from(a: Dict[str, float], flow_per_day: float, capex: float = None, include_demand: bool = True) -> DwptLcos: """前提辞書から DWPT 総 LCOS を計算する共通入口（受電器按分は a['dwpt_km_per_vehicle_year']）。""" e1 = energy_per_pass_kwh(a["peak_power_kw"], 1000.0, a["speed_kmh"], a["coupling_duty"]) return dwpt_lcos( a["capex_dwpt_yen_per_lane_km"] if capex is None else capex, a["discount_rate"], a["lifetime_years"], a["om_fraction"], flow_per_day, a["peak_power_kw"], a["speed_kmh"], a["coupling_duty"], a["grid_to_battery_eff"], a["electricity_yen_per_kwh"], a["receiver_cost_yen"], a["receiver_lifetime_years"], e1 * a.get("dwpt_km_per_vehicle_year", DWPT_KM_PER_VEHICLE_YEAR_CORRIDOR), a["demand_charge_yen_per_kw_month"] if include_demand else 0.0, a.get("dwpt_contract_kw_min", 150.0), a.get("peak_hour_share", 0.10), ) # ============================================================================= # 3. 代替案 A：電池を積む（電池スループット費用） # ============================================================================= @dataclass class BatteryCase: delta_kwh: float capex_yen: float mass_kg: float cycles_per_year: float throughput_yen_per_kwh: float payload_penalty_yen_per_kwh: float energy_penalty_yen_per_kwh: float total_yen_per_kwh: float def battery_throughput_cost( delta_kwh: float, pack_yen_per_kwh: float, cycle_life: float, truck_life_years: float, rate: float, truck_kwh_per_km: float, truck_km_per_year: float, base_pack_kwh: float, kg_per_kwh: float, payload_value_yen_per_tkm: float, payload_binding_fraction: float, electricity_yen_per_kwh: float, energy_penalty_per_tonne: float = 0.012, delta_utilisation: float = 1.0, ) -> BatteryCase: """「DWPT で道路から受け取る代わりに、電池を ΔE kWh 余分に積んで静止充電で賄う」費用を 1 kWh スループットあたり [円/kWh] で表す。 - 電池年額 = ΔE × 単価 × CRF（車両寿命かサイクル寿命の短い方で償却） - 年間スループット = ΔE × 年間等価フルサイクル数（年間使用 kWh ÷ パック容量で上限） - 質量ペナルティ = Δm[t] × 運賃[円/t-km] × 積載が重量制約で縛られる割合 ÷ 電費[kWh/km] - 電費ペナルティ = Δm[t] × (1 t あたり電費増 ≈ 1.2 %) × 電費 × 電気代 ÷ 電費""" annual_kwh_use = truck_kwh_per_km * truck_km_per_year pack_kwh = base_pack_kwh + delta_kwh cycles_per_year = annual_kwh_use / pack_kwh life_by_cycles = cycle_life / cycles_per_year life = min(truck_life_years, life_by_cycles) capex = delta_kwh * pack_yen_per_kwh annual_cost = capex * crf(rate, life) throughput = delta_kwh * cycles_per_year * delta_utilisation # ΔE が担う年間 kWh（限界 kWh の実使用率を掛ける） thr = annual_cost / throughput if throughput > 0 else float("inf") mass_kg = delta_kwh * kg_per_kwh # 質量ペナルティは全走行 km で発生する年額。分母は ΔE が担う年間スループット（償却と同じ分母）に揃える。 # 旧版は車両の全消費 kWh で割っており、ΔE/パック容量 倍だけ過小だった（レビュー R3 指摘で修正）。 payload_pen_annual = (mass_kg / 1000.0) * payload_value_yen_per_tkm * payload_binding_fraction * truck_km_per_year payload_pen = payload_pen_annual / throughput if throughput > 0 else float("inf") energy_pen_annual = (mass_kg / 1000.0) * energy_penalty_per_tonne * annual_kwh_use * electricity_yen_per_kwh energy_pen = energy_pen_annual / throughput if throughput > 0 else float("inf") return BatteryCase(delta_kwh, capex, mass_kg, cycles_per_year, thr, payload_pen, energy_pen, thr + payload_pen + energy_pen) # ============================================================================= # 4. 代替案 B：MCS（メガワット充電） # ============================================================================= @dataclass class McsLcos: infra_yen_per_kwh: float demand_charge_yen_per_kwh: float electricity_yen_per_kwh: float time_yen_per_kwh: float total_yen_per_kwh: float annual_kwh_per_stall: float def mcs_lcos( charger_capex_yen: float, grid_capex_yen_per_stall: float, rate: float, lifetime_years: float, om_fraction: float, utilisation: float, power_kw: float, avg_power_fraction: float, electricity_yen_per_kwh: float, demand_charge_yen_per_kw_month: float, eff: float, time_cost_yen_per_h: float, logistics_overlap: float, coincidence: float = 0.6, ) -> McsLcos: """MCS 1 口あたり LCOS [円/kWh]。 - 年間供給 = P × 8,760 h × 時間稼働率 × 充電中平均出力比（SOC テーパ） - infra = (充電器＋系統連系の年額) ÷ 年間供給 - 基本料金 = 契約 kW（= P × 同時率）× 月額 × 12 ÷ 年間供給 - 時間費用 = (運転手＋車両の時間価値 ÷ 平均充電出力) × (1 − 休憩等との重なり率)""" annual_kwh = power_kw * HOURS_PER_YEAR * utilisation * avg_power_fraction ann = annualised_cost(charger_capex_yen + grid_capex_yen_per_stall, rate, lifetime_years, om_fraction) infra = ann / annual_kwh demand = demand_charge_yen_per_kw_month * 12.0 * power_kw * coincidence / annual_kwh elec = electricity_yen_per_kwh / eff time_cost = time_cost_yen_per_h / (power_kw * avg_power_fraction) * (1.0 - logistics_overlap) return McsLcos(infra, demand, elec, time_cost, infra + demand + elec + time_cost, annual_kwh) def mcs_charge_time_min( energy_kwh: float, power_kw: float, pack_kwh: float, soc_start: float = 0.15, max_c_rate: float = 2.0, taper_start: float = 0.70, taper_end_fraction: float = 0.25, dt_s: float = 1.0, ) -> float: """SOC カーブ付きの充電時間 [min]。 出力 = min(P_charger, C レート上限 × 容量)。SOC 70 % までは一定、以降 100 % で 25 % へ線形テーパ。""" p_cap = min(power_kw, max_c_rate * pack_kwh) soc = soc_start soc_target = min(1.0, soc_start + energy_kwh / pack_kwh) t = 0.0 while soc < soc_target and t < 36000: if soc < taper_start: p = p_cap else: frac = (soc - taper_start) / (1.0 - taper_start) p = p_cap * (1.0 - frac * (1.0 - taper_end_fraction)) soc += p * dt_s / SEC_PER_HOUR / pack_kwh t += dt_s return t / 60.0 # ============================================================================= # 5. MCS ハブの電力集中（DWPT 側の有力反論） # ============================================================================= @dataclass class HubCase: simultaneous_trucks: int peak_mw: float grid_connection_yen: float buffer_battery_kwh: float buffer_battery_yen: float total_yen: float yen_per_kwh_delivered: float def mcs_hub( simultaneous_trucks: int, power_kw: float, grid_yen_per_kw: float, buffer_hours: float, pack_yen_per_kwh: float, daily_sessions_per_stall: float, kwh_per_session: float, rate: float, life_years: float, om_fraction: float, ) -> HubCase: """N 台同時 MCS 充電ハブの系統連系費と、ピークを均す電池バッファ費。 ピーク MW = N × P。連系費 = ピーク kW × 円/kW。バッファ = ピーク × 時間 × 単価。""" peak_kw = simultaneous_trucks * power_kw grid = peak_kw * grid_yen_per_kw buf_kwh = peak_kw * buffer_hours buf_yen = buf_kwh * pack_yen_per_kwh total = grid + buf_yen annual_kwh = simultaneous_trucks * daily_sessions_per_stall * kwh_per_session * DAYS_PER_YEAR ann = annualised_cost(total, rate, life_years, om_fraction) return HubCase(simultaneous_trucks, peak_kw / 1000.0, grid, buf_kwh, buf_yen, total, ann / annual_kwh) # ============================================================================= # 6. ケース計算 # ============================================================================= @dataclass class TateyamaCase: year: int ev_fraction: float dwpt_penetration: float target_flow_per_day: float e_pass_kwh: float pass_time_s: float coil_seconds_per_day: float coil_minutes_per_day: float coil_hours_per_year: float coil_time_share: float annual_kwh_300m: float kwh_per_m_year: float mwh_per_km_year: float capacity_utilisation: float lcos_infra_yen_per_kwh: float lcos_total_yen_per_kwh: float lcos_receiver_if_only_this_300m_yen_per_kwh: float lcos_demand_charge_yen_per_kwh: float = 0.0 lcos_total_excl_demand_yen_per_kwh: float = 0.0 def tateyama_grid(a: Dict[str, float], years=(2026, 2030, 2035, 2040), pens=(0.01, 0.05, 0.10, 0.25, 0.50, 1.0), dwpt_km_per_vehicle_year: float = DWPT_KM_PER_VEHICLE_YEAR_CORRIDOR) -> List[TateyamaCase]: """館山道 300 m 区間の年×搭載率グリッド。 受電器費用は「対象車が年間 dwpt_km_per_vehicle_year km の DWPT 区間を走る」前提で kWh に割り付ける （既定 25,000 km/年＝幹線回廊が整備済みという DWPT 有利側の仮定）。『この 300 m しか存在しない』場合の受電器 LCOS は別列に出す。""" out: List[TateyamaCase] = [] D = a["section_length_m"] for y in years: evf = a[f"ev_fraction_heavy_{y}"] for pen in pens: flow = target_flow_per_day(a["aadt_total"], a["heavy_vehicle_fraction"], evf, pen, a["lane_fraction"]) e_pass = energy_per_pass_kwh(a["peak_power_kw"], D, a["speed_kmh"], a["coupling_duty"]) t_pass = pass_time_s(D, a["speed_kmh"]) occ = occupancy_time_s(D, a["vehicle_length_m"], a["speed_kmh"]) sec_day = flow * occ annual = annual_energy_kwh(flow, e_pass) l = dwpt_lcos_from(dict(a, dwpt_km_per_vehicle_year=dwpt_km_per_vehicle_year), flow) recv_only_here = annualised_cost(a["receiver_cost_yen"], a["discount_rate"], a["receiver_lifetime_years"], 0.0) / (e_pass * DAYS_PER_YEAR) out.append(TateyamaCase( y, evf, pen, flow, e_pass, t_pass, sec_day, sec_day / 60.0, sec_day * DAYS_PER_YEAR / 3600.0, coil_time_share(flow, D, a["vehicle_length_m"], a["speed_kmh"]), annual, kwh_per_m_year(annual, D), mwh_per_km_year(annual, D), power_capacity_utilisation(annual, a["peak_power_kw"]), l.infra_yen_per_kwh, l.total_yen_per_kwh, recv_only_here, l.demand_charge_yen_per_kwh, l.total_yen_per_kwh - l.demand_charge_yen_per_kwh, )) return out def break_even_capex(a: Dict[str, float], flow_per_day: float, alternative_yen_per_kwh: float, dwpt_km_per_vehicle_year: float = DWPT_KM_PER_VEHICLE_YEAR_CORRIDOR) -> float: """与えた対象車流量で DWPT 総 LCOS が代替案と等しくなる CAPEX [円/lane-km]。 インフラ LCOS は CAPEX に比例するので閉形式で解ける。""" e1 = energy_per_pass_kwh(a["peak_power_kw"], 1000.0, a["speed_kmh"], a["coupling_duty"]) annual_kwh = annual_energy_kwh(flow_per_day, e1) recv = annualised_cost(a["receiver_cost_yen"], a["discount_rate"], a["receiver_lifetime_years"], 0.0) / (e1 * a.get("dwpt_km_per_vehicle_year", dwpt_km_per_vehicle_year)) elec = a["electricity_yen_per_kwh"] / a["grid_to_battery_eff"] demand = dwpt_lcos_from(a, flow_per_day).demand_charge_yen_per_kwh room = alternative_yen_per_kwh - recv - elec - demand if room <= 0: return 0.0 return room * annual_kwh / (crf(a["discount_rate"], a["lifetime_years"]) + a["om_fraction"]) def kw_to_kwh_table(power_kw: float, v_kmh: float, lengths_m=(300, 1000, 5000, 10000, 100000)) -> List[dict]: rows = [] for D in lengths_m: rows.append({ "length_m": D, "speed_kmh": v_kmh, "power_kw": power_kw, "pass_time_s": round(pass_time_s(D, v_kmh), 1), "kwh_per_pass_duty100": round(energy_per_pass_kwh(power_kw, D, v_kmh, 1.0), 3), "kwh_per_pass_duty75": round(energy_per_pass_kwh(power_kw, D, v_kmh, 0.75), 3), "truck_km_equivalent_at_1p3kwh_per_km": round(energy_per_pass_kwh(power_kw, D, v_kmh, 0.75) / 1.3, 2), }) return rows def speed_sweep(power_kw: float, section_m: float = 100.0, speeds=(120, 100, 80, 50, 20, 10, 5), stop_seconds=(0, 60, 300, 1800)) -> List[dict]: """同じ 100 m で速度を落とすと受電量がどう増えるか。最後に「停止」を加える。""" rows = [] for v in speeds: t = pass_time_s(section_m, v) rows.append({"mode": f"{v} km/h", "speed_kmh": v, "time_on_coil_s": round(t, 1), "kwh": round(power_kw * t / 3600, 3)}) for s in stop_seconds[1:]: rows.append({"mode": f"stop {s}s", "speed_kmh": 0, "time_on_coil_s": s, "kwh": round(power_kw * s / 3600, 3)}) return rows def length_scaling(a: Dict[str, float], flow_per_day: float, lengths_m=(300, 1000, 3000, 10000, 30000)) -> List[dict]: """区間長 n 倍 → kWh/pass n 倍、CAPEX n 倍なら kWh/m/year と LCOS は不変（比例モデル）。 固定費（受変電・制御・通信）と車両の受電上限を入れた列も併記する：固定費は短区間を不利に、受電上限は長区間を不利にする。""" rows = [] fixed = a.get("fixed_capex_per_segment_yen", 0.0) usable = a.get("usable_kwh_per_km", 1e9) for D in lengths_m: e_pass = energy_per_pass_kwh(a["peak_power_kw"], D, a["speed_kmh"], a["coupling_duty"]) annual = annual_energy_kwh(flow_per_day, e_pass) capex = a["capex_dwpt_yen_per_lane_km"] * D / 1000.0 ann = annualised_cost(capex, a["discount_rate"], a["lifetime_years"], a["om_fraction"]) e_cap = min(e_pass, usable * D / 1000.0) annual_cap = annual_energy_kwh(flow_per_day, e_cap) ann_fixed = annualised_cost(capex + fixed, a["discount_rate"], a["lifetime_years"], a["om_fraction"]) rows.append({"length_m": D, "kwh_per_pass": round(e_pass, 2), "annual_kwh": round(annual), "capex_yen": round(capex), "annualised_yen": round(ann), "kwh_per_m_year": round(annual / D, 2), "infra_yen_per_kwh": round(ann / annual, 1) if annual else None, "infra_yen_per_kwh_with_fixed_cost_and_cap": round(ann_fixed / annual_cap, 1) if annual_cap else None, "effective_capex_yen_per_lane_km": round((capex + fixed) / (D / 1000.0))}) return rows def break_even_flow(a: Dict[str, float], alternative_yen_per_kwh: float) -> float: """DWPT の総 LCOS が代替案と等しくなる対象車流量 [台/日・敷設車線] を二分法で求める。""" def total(flow: float) -> float: return dwpt_lcos_from(a, flow).total_yen_per_kwh lo, hi = 1.0, 1e7 if total(hi) > alternative_yen_per_kwh: return float("inf") for _ in range(200): mid = math.sqrt(lo * hi) if total(mid) > alternative_yen_per_kwh: lo = mid else: hi = mid return hi def alternative_lcos(a: Dict[str, float], year: int = 2026, delta_kwh: float = 200.0, base_pack_kwh: float = 400.0, payload_binding: float = 0.3, battery_price: float = None) -> dict: """「電池 ΔE を積む ＋ MCS で静止充電」の合計 LCOS を返す。""" b = battery_throughput_cost( delta_kwh, a[f"battery_pack_yen_per_kwh_{year}"] if battery_price is None else battery_price, a["battery_cycle_life"], a["truck_life_years"], a["discount_rate"], a["truck_kwh_per_km"], a["truck_km_per_year"], base_pack_kwh, a["battery_kg_per_kwh"], a["payload_value_yen_per_tonne_km"], payload_binding, a["electricity_yen_per_kwh"], delta_utilisation=a.get("delta_kwh_utilisation", 1.0), ) m = mcs_lcos( a["mcs_charger_capex_yen"], a["mcs_site_grid_capex_yen"] / a["mcs_stalls_per_site"], a["discount_rate"], a["mcs_lifetime_years"], a["mcs_om_fraction"], a["mcs_utilization"], a["mcs_power_kw"], 0.7, a["electricity_yen_per_kwh"], a["demand_charge_yen_per_kw_month"], 0.93, a["driver_vehicle_hour_yen"], a["logistics_overlap"], ) return {"battery": b, "mcs": m, "total_yen_per_kwh": b.total_yen_per_kwh + m.total_yen_per_kwh} def depot_charging_lcos(a: Dict[str, float], year: int = 2026, delta_kwh: float = 200.0) -> dict: """短距離運行の代替案：車庫での夜間普通充電（時間費用ゼロ、基本料金込み単価）＋ 電池 ΔE。""" b = battery_throughput_cost( delta_kwh, a[f"battery_pack_yen_per_kwh_{year}"], a["battery_cycle_life"], a["truck_life_years"], a["discount_rate"], a["truck_kwh_per_km"], a["truck_km_per_year"], 400.0, a["battery_kg_per_kwh"], a["payload_value_yen_per_tonne_km"], 0.3, a["electricity_yen_per_kwh"], delta_utilisation=a.get("delta_kwh_utilisation", 1.0), ) charger_kw = 100.0 annual_kwh = charger_kw * HOURS_PER_YEAR * 0.20 # 夜間 8 h × 60 % 稼働 ≈ 20 % infra = annualised_cost(a.get("depot_charger_capex_yen", 5e6), a["discount_rate"], 12.0, 0.03) / annual_kwh elec = a.get("depot_electricity_yen_per_kwh", 28.0) / 0.93 return {"battery": b, "infra_yen_per_kwh": infra, "electricity_yen_per_kwh": elec, "total_yen_per_kwh": b.total_yen_per_kwh + infra + elec} # ============================================================================= # 6b. 幹線ケース：東名・新東名・東北道で DWPT が並ぶのに必要な受電器搭載率 # ============================================================================= TRUNK_LINES = [ # (名称, 大型車 24h 上下計 [台/日], 車線数(上下計), 出典) ("館山道 君津–君津PAスマート", 1699, 4, "R3 センサス"), ("東名 厚木–秦野中井", 49469, 6, "R3 センサス"), ("新東名 御殿場JCT–長泉沼津", 35128, 6, "R3 センサス"), ("東北道 久喜–加須", 32436, 6, "R3 センサス"), ] def trunk_line_cases(a: Dict[str, float], low: Dict[str, float], high: Dict[str, float], lane_share_of_heavy: float = 0.9) -> List[dict]: """幹線で『敷設車線 1 本』に流れる大型 BEV 受電車の流量と、損益分岐に必要な搭載率。 大型車は方向別に半分、そのうち走行車線（敷設車線）を使う割合 lane_share_of_heavy。""" rows = [] for name, heavy_both, lanes, src in TRUNK_LINES: heavy_dir_lane = heavy_both / 2.0 * lane_share_of_heavy for capex_name, capex in (("LOW", low["capex_dwpt_yen_per_lane_km"]), ("BASE", a["capex_dwpt_yen_per_lane_km"])): aa = dict(a, capex_dwpt_yen_per_lane_km=capex) alt = alternative_lcos(aa, 2040)["total_yen_per_kwh"] be = break_even_flow(aa, alt) for ev_name, ev in (("BASE 20%", a["ev_fraction_heavy_2040"]), ("HIGH 45%", high["ev_fraction_heavy_2040"])): bev_lane = heavy_dir_lane * ev req_pen = be / bev_lane if bev_lane > 0 else float("inf") l30 = dwpt_lcos_from(aa, bev_lane * 0.3) rows.append({"road": name, "heavy_both_dir_per_day": heavy_both, "heavy_per_lane_dir": round(heavy_dir_lane), "capex_case": capex_name, "bev_2040_case": ev_name, "bev_per_lane_dir_2040": round(bev_lane), "alt_yen_per_kwh": round(alt, 1), "break_even_flow_per_lane": round(be), "required_penetration": round(req_pen, 2) if req_pen != float("inf") else "inf", "dwpt_total_yen_per_kwh_at_pen30": round(l30.total_yen_per_kwh, 1), "source": src}) return rows def financing_factor(rate: float, life: float, om: float) -> float: return crf(rate, life) + om def split_financing_factor(civil_share: float = 0.5, civil_rate: float = 0.04, civil_life: float = 50, elec_rate: float = 0.04, elec_life: float = 15, om: float = 0.02) -> float: """土木部分（配管・筐体・舗装構造）を道路本体と同じ 50 年・4%、電子機器を 15 年・4% で分けて償却した年額係数。""" return civil_share * crf(civil_rate, civil_life) + (1 - civil_share) * crf(elec_rate, elec_life) + om # ============================================================================= # 7. Dream Case（DWPT に極端に有利な前提） # ============================================================================= def dream_case_assumptions(low: Dict[str, float], base: Dict[str, float], high: Dict[str, float]) -> Dict[str, float]: d = dict(base) d.update({ "capex_dwpt_yen_per_lane_km": low["capex_dwpt_yen_per_lane_km"], "lifetime_years": low["lifetime_years"], "om_fraction": low["om_fraction"], "discount_rate": low["discount_rate"], "dwpt_penetration": 1.0, "aadt_total": high["aadt_total"], "heavy_vehicle_fraction": high["heavy_vehicle_fraction"], "ev_fraction_heavy_2040": high["ev_fraction_heavy_2040"], "coupling_duty": high["coupling_duty"], "grid_to_battery_eff": high["grid_to_battery_eff"], "speed_kmh": high["speed_kmh"], # 80 km/h（遅いほど有利） "lane_fraction": high["lane_fraction"], "receiver_cost_yen": low["receiver_cost_yen"], "receiver_lifetime_years": low["receiver_lifetime_years"], "peak_power_kw": 300.0, # 受電器多重化（VINCI A10 でピーク 300 kW 超の公道実測。レビュー R1-01/R4-17） "dwpt_contract_kw_min": low["dwpt_contract_kw_min"], "fixed_capex_per_segment_yen": low["fixed_capex_per_segment_yen"], # 代替案は不利側 "battery_pack_yen_per_kwh_2026": high["battery_pack_yen_per_kwh_2026"], "battery_pack_yen_per_kwh_2040": high["battery_pack_yen_per_kwh_2040"], "battery_cycle_life": low["battery_cycle_life"], "battery_kg_per_kwh": high["battery_kg_per_kwh"], "mcs_charger_capex_yen": high["mcs_charger_capex_yen"], "mcs_site_grid_capex_yen": high["mcs_site_grid_capex_yen"], "mcs_utilization": low["mcs_utilization"], "logistics_overlap": 0.0, "driver_vehicle_hour_yen": high["driver_vehicle_hour_yen"], "payload_value_yen_per_tonne_km": high["payload_value_yen_per_tonne_km"], }) return d # ============================================================================= # 8. 実行 # ============================================================================= def write_csv(path: pathlib.Path, rows: Iterable[dict]) -> None: rows = list(rows) if not rows: return path.parent.mkdir(parents=True, exist_ok=True) with open(path, "w", newline="", encoding="utf-8") as f: w = csv.DictWriter(f, fieldnames=list(rows[0].keys())) w.writeheader() w.writerows(rows) def main() -> None: low, base, high = load_assumptions("low"), load_assumptions("base"), load_assumptions("high") RESULTS.mkdir(exist_ok=True) # (1) kW → kWh write_csv(RESULTS / "kw_to_kwh.csv", kw_to_kwh_table(base["peak_power_kw"], base["speed_kmh"])) # (2) 館山道グリッド（BASE） grid = tateyama_grid(base) write_csv(RESULTS / "tateyama_grid_base.csv", [asdict(c) for c in grid]) write_csv(RESULTS / "tateyama_grid_high_traffic.csv", [asdict(c) for c in tateyama_grid(dict(base, aadt_total=high["aadt_total"], heavy_vehicle_fraction=high["heavy_vehicle_fraction"], lane_fraction=high["lane_fraction"]))]) # (3) 速度掃引・長さスケーリング write_csv(RESULTS / "speed_sweep.csv", speed_sweep(base["peak_power_kw"])) flow_2040_full = target_flow_per_day(base["aadt_total"], base["heavy_vehicle_fraction"], base["ev_fraction_heavy_2040"], 1.0, base["lane_fraction"]) write_csv(RESULTS / "length_scaling.csv", length_scaling(base, flow_2040_full)) # (4) 電池ケース rows = [] for y in (2026, 2030, 2035, 2040): for dk in (50, 100, 200, 300, 500, 800): b = battery_throughput_cost(dk, base[f"battery_pack_yen_per_kwh_{y}"], base["battery_cycle_life"], base["truck_life_years"], base["discount_rate"], base["truck_kwh_per_km"], base["truck_km_per_year"], 400.0, base["battery_kg_per_kwh"], base["payload_value_yen_per_tonne_km"], 0.3, base["electricity_yen_per_kwh"]) rows.append({"year": y, **{k: (round(v, 3) if isinstance(v, float) else v) for k, v in asdict(b).items()}}) write_csv(RESULTS / "battery_cases.csv", rows) # (5) MCS 充電時間 rows = [] for p in (1000, 1500, 3000): for pack in (500, 600, 900): rows.append({"mcs_kw": p, "pack_kwh": pack, "minutes_for_500kwh": round(mcs_charge_time_min(500, p, pack), 1), "minutes_for_300kwh": round(mcs_charge_time_min(300, p, pack), 1)}) write_csv(RESULTS / "mcs_charge_time.csv", rows) # (6) MCS LCOS と overlap 掃引 rows = [] for ov in (0.0, 0.25, 0.5, 0.75, 1.0): for util in (0.05, 0.15, 0.35): m = mcs_lcos(base["mcs_charger_capex_yen"], base["mcs_site_grid_capex_yen"] / base["mcs_stalls_per_site"], base["discount_rate"], base["mcs_lifetime_years"], base["mcs_om_fraction"], util, base["mcs_power_kw"], 0.7, base["electricity_yen_per_kwh"], base["demand_charge_yen_per_kw_month"], 0.93, base["driver_vehicle_hour_yen"], ov) rows.append({"overlap": ov, "utilisation": util, **{k: round(v, 2) for k, v in asdict(m).items()}}) write_csv(RESULTS / "mcs_lcos.csv", rows) # (7) MCS ハブ電力集中 rows = [] for n in (10, 20, 50, 100): h = mcs_hub(n, base["mcs_power_kw"], 30000.0, 1.0, base["battery_pack_yen_per_kwh_2026"], 8.0, 400.0, base["discount_rate"], 15.0, 0.02) rows.append({k: (round(v, 2) if isinstance(v, float) else v) for k, v in asdict(h).items()}) write_csv(RESULTS / "mcs_hub.csv", rows) # (8) 損益分岐流量 rows = [] for case_name, a in (("low_capex", dict(base, capex_dwpt_yen_per_lane_km=low["capex_dwpt_yen_per_lane_km"])), ("base", base), ("high_capex", dict(base, capex_dwpt_yen_per_lane_km=high["capex_dwpt_yen_per_lane_km"]))): for y in (2026, 2030, 2040): alt = alternative_lcos(a, y) be = break_even_flow(a, alt["total_yen_per_kwh"]) be_nd = break_even_flow(dict(a, demand_charge_yen_per_kw_month=0.0), alt["total_yen_per_kwh"]) rows.append({"case": case_name, "year": y, "alternative_yen_per_kwh": round(alt["total_yen_per_kwh"], 1), "battery_component": round(alt["battery"].total_yen_per_kwh, 1), "mcs_component": round(alt["mcs"].total_yen_per_kwh, 1), "break_even_target_vehicles_per_day_per_lane": round(be) if be != float("inf") else "inf", "break_even_excl_demand_charge": round(be_nd) if be_nd != float("inf") else "inf", "capex_yen_per_lane_km": a["capex_dwpt_yen_per_lane_km"]}) write_csv(RESULTS / "break_even.csv", rows) # (7b) 車庫充電代替 rows = [] for y in (2026, 2030, 2040): dp = depot_charging_lcos(base, y) rows.append({"year": y, "battery_yen_per_kwh": round(dp["battery"].total_yen_per_kwh, 1), "depot_infra_yen_per_kwh": round(dp["infra_yen_per_kwh"], 1), "depot_electricity_yen_per_kwh": round(dp["electricity_yen_per_kwh"], 1), "total_yen_per_kwh": round(dp["total_yen_per_kwh"], 1), "mcs_alternative_total_yen_per_kwh": round(alternative_lcos(base, y)["total_yen_per_kwh"], 1)}) write_csv(RESULTS / "depot_alternative.csv", rows) # (8a) 幹線ケース write_csv(RESULTS / "trunk_line_cases.csv", trunk_line_cases(base, low, high)) # (8a') 評価期間・割引率の感度：年額係数 rows = [] for name, fac in (("BASE 20y 5% O&M2%", financing_factor(0.05, 20, 0.02)), ("30y 3% O&M1%", financing_factor(0.03, 30, 0.01)), ("道路CBA 50y 4% O&M2%（全額）", financing_factor(0.04, 50, 0.02)), ("土木50%@50y4% + 機器50%@15y4% O&M2%", split_financing_factor()), ("10y 8% O&M5%", financing_factor(0.08, 10, 0.05))): flow = target_flow_per_day(base["aadt_total"], base["heavy_vehicle_fraction"], base["ev_fraction_heavy_2040"], 1.0, base["lane_fraction"]) e_pass = energy_per_pass_kwh(base["peak_power_kw"], base["section_length_m"], base["speed_kmh"], base["coupling_duty"]) annual = annual_energy_kwh(flow, e_pass) ann = base["capex_dwpt_yen_per_lane_km"] * base["section_length_m"] / 1000.0 * fac rows.append({"financing": name, "annual_factor": round(fac, 4), "annualised_300m_yen": round(ann), "infra_lcos_2040_pen100_yen_per_kwh": round(ann / annual, 1), "relative_to_base": round(fac / financing_factor(0.05, 20, 0.02), 3)}) write_csv(RESULTS / "financing_sensitivity.csv", rows) # (8b) 損益分岐 CAPEX：館山道の各シナリオ流量で、いくらまで CAPEX が下がれば並ぶか rows = [] for fin_name, fin in (("base_20y_5pct_om2", {}), ("best_30y_3pct_om1", {"lifetime_years": low["lifetime_years"], "discount_rate": low["discount_rate"], "om_fraction": low["om_fraction"]})): a = dict(base, **fin) for y in (2030, 2035, 2040): for pen in (0.1, 0.5, 1.0): flow = target_flow_per_day(a["aadt_total"], a["heavy_vehicle_fraction"], a[f"ev_fraction_heavy_{y}"], pen, a["lane_fraction"]) alt = alternative_lcos(a, y)["total_yen_per_kwh"] rows.append({"financing": fin_name, "year": y, "dwpt_penetration": pen, "target_flow_per_day": round(flow, 1), "alternative_yen_per_kwh": round(alt, 1), "break_even_capex_yen_per_lane_km": round(break_even_capex(a, flow, alt)), "break_even_capex_oku_yen": round(break_even_capex(a, flow, alt) / 1e8, 2)}) write_csv(RESULTS / "break_even_capex.csv", rows) # (9) Dream Case d = dream_case_assumptions(low, base, high) flow_dream = target_flow_per_day(d["aadt_total"], d["heavy_vehicle_fraction"], d["ev_fraction_heavy_2040"], 1.0, d["lane_fraction"]) l = dwpt_lcos_from(d, flow_dream) alt = alternative_lcos(d, 2040) dream = {"target_flow_per_day": round(flow_dream), "dwpt_infra_yen_per_kwh": round(l.infra_yen_per_kwh, 1), "dwpt_total_yen_per_kwh": round(l.total_yen_per_kwh, 1), "alternative_total_yen_per_kwh": round(alt["total_yen_per_kwh"], 1), "battery_component": round(alt["battery"].total_yen_per_kwh, 1), "mcs_component": round(alt["mcs"].total_yen_per_kwh, 1), "annual_mwh_per_lane_km": round(l.annual_kwh_per_lane_km / 1000, 1), "break_even_flow": break_even_flow(d, alt["total_yen_per_kwh"])} with open(RESULTS / "dream_case.json", "w", encoding="utf-8") as f: json.dump({"assumptions": d, "result": dream}, f, ensure_ascii=False, indent=2) print(json.dumps(dream, ensure_ascii=False, indent=2)) for c in grid: if c.year in (2026, 2040) and c.dwpt_penetration in (0.1, 1.0): print(f"{c.year} pen={c.dwpt_penetration:.0%}: flow={c.target_flow_per_day:.1f}/day, coil={c.coil_minutes_per_day:.1f} min/day, " f"kWh/m/yr={c.kwh_per_m_year:.2f}, LCOS infra={c.lcos_infra_yen_per_kwh:,.0f} 円/kWh") if __name__ == "__main__": main()