#!/usr/bin/env python3 """Delivered-cost model, version 2.1 (October 24, 2026). Compares three power-supply arrangements for a 200 MW IT campus on a common reliability baseline, with the grid bill built from Dominion Energy Virginia tariff components (see TARIFF dict; sources in the paper's Section 12 and the workbook's Tariff sheet). Arrangements A Grid + diesel backup (N+1), common baseline C Grid + dual-fuel gas engines (N+2) + grid-forming battery, common baseline E Onsite engines only (islanded), common baseline Common baseline (identical in all three, so it cancels in differences but is included so absolute $/MWh are like-for-like): 5-minute lithium-ion UPS on the IT and CDU load; ride-through-rated cooling drives, tuned transfer settings and the electromagnetic-transient study; 48-72 h liquid-fuel storage for whatever backup fleet is installed. All figures are 2026 US dollars, unsubsidized, engineering estimates unless a tariff component. Capital is annualized at a 10% capital recovery factor. """ import json, sys IT_MW = 200.0 PUE = 1.15 AVG_MW = IT_MW * PUE # 230 MW average facility load PEAK_MW = 240.0 # billing demand, kW basis MWH = AVG_MW * 8760.0 # 2,014,800 MWh/yr CRF = 0.10 # ---------------- Tariff inputs (filled from research; see TARIFF_NOTES) ---- # Units: demand charges $/kW-month; energy and riders $/kWh; customer $/month. TARIFF = dict( customer_charge_month=0.0, dist_demand_kw_month=0.0, # on T&D billing demand (85% minimum) trans_demand_kw_month=0.0, # on T&D billing demand gen_demand_kw_month=0.0, # on generation billing demand (60% minimum) riders_trans_kw_month=0.0, # Rider T1, on T&D billing demand riders_gen_kw_month=0.0, # Riders E, GEN, SNA, CERC, SMR, CE, OSW, on generation billing demand energy_kwh=0.0, # base generation + distribution energy charges fuel_rider_kwh=0.0, # Rider A other_riders_kwh=0.0, # CCR, RPS, C1A, C4A, RGGI, PIPP, deferred fuel offpeak_gen_demand_kw_month=0.0, # GS-4 II.B.2, on off-peak demand above 90% of on-peak demand (Section VII) offpeak_threshold=0.90, billing_days=365, # GS-4 Section IX prorates the schedule's charges by days/30 cap_share_of_gen=1.0, # share of the generation demand charge scaled in the uncapped case fuel_mult_high=1.5, # Rider A multiplier in the high-gas case billing_demand_mw=PEAK_MW, # fallback: if components are zero, use these flat assumptions (version 1) flat_energy_mwh=65.0, flat_capacity_mw_yr=121705.0, ) def grid_bill(t, mwh=MWH, demand_mw=PEAK_MW, demand_credit_mw=0.0): """Annual retail bill in $ from tariff components. If component rates are all zero, fall back to the version-1 flat energy + PJM capacity assumption.""" comp = (t['dist_demand_kw_month'] + t['trans_demand_kw_month'] + t['gen_demand_kw_month'] + t['energy_kwh'] + t['fuel_rider_kwh']) t = dict(t); t.setdefault('fuel_mult', 1.0); t.setdefault('gen_mult', 1.0) if comp == 0: return t['flat_energy_mwh'] * mwh + (demand_mw - demand_credit_mw) * t['flat_capacity_mw_yr'] kw = demand_mw * 1000.0 kw_gen = max(demand_mw - demand_credit_mw, 0.60 * demand_mw) * 1000.0 # GS-5 60% generation minimum kw_td = max(demand_mw, 0.85 * demand_mw) * 1000.0 # 85% T&D minimum kw_off = max(kw - t['offpeak_threshold'] * kw_gen, 0.0) # off-peak demand above 90% of on-peak kwh = mwh * 1000.0 m_sched = t['billing_days'] / 30.0 # Section IX proration of GS-4 charges return (t['customer_charge_month'] * m_sched + (t['dist_demand_kw_month'] + t['trans_demand_kw_month']) * kw_td * m_sched + t['riders_trans_kw_month'] * kw_td * 12 + t['gen_demand_kw_month'] * t['gen_mult'] * kw_gen * m_sched + t['offpeak_gen_demand_kw_month'] * kw_off * m_sched + t['riders_gen_kw_month'] * kw_gen * 12 + (t['energy_kwh'] + t['fuel_rider_kwh'] * t['fuel_mult'] + t['other_riders_kwh']) * kwh) # ---------------- Equipment and fuel assumptions -------------------------- P = dict( ups_capex_kw=450.0, ups_om_kw_yr=5.0, ups_mw=210.0, # IT + CDU pumps ridethrough_pkg_capex=15e6, # drives, controls, EMT study (one-off) fuel_tank_capex_per_mw_backup=40e3, # 72 h liquid fuel storage, $/MW of backup diesel_mw=264.0, diesel_capex_kw=800.0, diesel_om_kw_yr=8.0, diesel_test_h=60.0, diesel_fuel_mwh=200.0, # test/emergency hours, $/MWh of diesel fuel eng_mw=320.0, eng_capex_kw=2000.0, eng_om_kw_yr=15.0, eng_heat_rate=8.5, eng_run_h=100.0, eng_vom_mwh=15.0, # outage/emergency hours per year gas_firm_transport_dth_day=0.40, # $/Dth-day reservation, firm gas_lateral_capex=25e6, bess_mw=120.0, bess_mwh=240.0, bess_capex_kwh=400.0, bess_om_kw_yr=10.0, bess_aug_frac=0.015, bess_cycles_yr=250, bess_rte=0.88, onsite_eng_mw=360.0, onsite_heat_rate=8.5, onsite_vom_mwh=15.0, ) CAP_CAPPED, CAP_UNCAPPED = 333.44, 542.83 # $/MW-day, PJM 2027/28 BRA clearing price and uncapped Dominion-zone simulation CASES = { 'Base': dict(gas=4.50, high_gas=False, uncapped=False), 'High gas ($8/MMBtu)': dict(gas=8.00, high_gas=True, uncapped=False), 'Uncapped capacity': dict(gas=4.50, high_gas=False, uncapped=True), } def baseline(): ups = P['ups_mw']*1e3*P['ups_capex_kw']*CRF + P['ups_mw']*1e3*P['ups_om_kw_yr'] rt = P['ridethrough_pkg_capex']*CRF return dict(ups=ups, ridethrough=rt) def tariff_for(case, t=TARIFF): tt=dict(t) tt['fuel_mult'] = t['fuel_mult_high'] if case['high_gas'] else 1.0 tt['gen_mult'] = 1.0 + t['cap_share_of_gen']*(CAP_UNCAPPED/CAP_CAPPED - 1.0) if case['uncapped'] else 1.0 if case['uncapped']: tt['flat_capacity_mw_yr'] = t['flat_capacity_mw_yr']*CAP_UNCAPPED*365/121705.0 if case['high_gas']: tt['flat_energy_mwh'] = t['flat_energy_mwh']*1.10 return tt def config_A(case, t=TARIFF): b = baseline() bill = grid_bill(tariff_for(case,t)) diesel = (P['diesel_mw']*1e3*P['diesel_capex_kw']*CRF + P['diesel_mw']*1e3*P['diesel_om_kw_yr'] + P['diesel_mw']*P['diesel_test_h']*P['diesel_fuel_mwh'] + P['diesel_mw']*P['fuel_tank_capex_per_mw_backup']*CRF) return dict(bill=bill, backup=diesel, **b) def config_C(case, t=TARIFF, credit_mw=0.0): b = baseline() bill = grid_bill(tariff_for(case,t), demand_credit_mw=credit_mw) eng = (P['eng_mw']*1e3*P['eng_capex_kw']*CRF + P['eng_mw']*1e3*P['eng_om_kw_yr'] + P['gas_lateral_capex']*CRF + P['eng_mw']*P['fuel_tank_capex_per_mw_backup']*CRF) dth_day = P['eng_mw']*P['eng_heat_rate']*24.0 # MMBtu/day at full load firm = dth_day*P['gas_firm_transport_dth_day']*365.0 run = P['eng_mw']*P['eng_run_h']*(P['eng_heat_rate']*case['gas'] + P['eng_vom_mwh']) bcap = P['bess_mwh']*1e3*P['bess_capex_kwh'] bess = bcap*CRF + P['bess_mw']*1e3*P['bess_om_kw_yr'] + bcap*P['bess_aug_frac'] tt=tariff_for(case,t) e_price = (tt['energy_kwh']+tt['fuel_rider_kwh']*tt['fuel_mult']+tt['other_riders_kwh'])*1000 if tt['energy_kwh'] else tt['flat_energy_mwh'] losses = P['bess_mwh']*P['bess_cycles_yr']*(1.0/P['bess_rte']-1.0)*e_price # charging energy above the energy delivered return dict(bill=bill, engines=eng, firm_gas=firm, engine_run=run, battery=bess+losses, **b) def config_E(case, t=TARIFF): b = baseline() eng = (P['onsite_eng_mw']*1e3*P['eng_capex_kw']*CRF + P['onsite_eng_mw']*1e3*P['eng_om_kw_yr'] + P['gas_lateral_capex']*CRF + P['onsite_eng_mw']*P['fuel_tank_capex_per_mw_backup']*CRF) fuel = MWH*(P['onsite_heat_rate']*case['gas'] + P['onsite_vom_mwh']) dth_day = PEAK_MW*P['onsite_heat_rate']*24.0 # firm transportation for the campus peak burn firm = dth_day*P['gas_firm_transport_dth_day']*365.0 bcap = P['bess_mwh']*1e3*P['bess_capex_kwh'] bess = bcap*CRF + P['bess_mw']*1e3*P['bess_om_kw_yr'] + bcap*P['bess_aug_frac'] losses = P['bess_mwh']*P['bess_cycles_yr']*(1.0/P['bess_rte']-1.0)*(P['onsite_heat_rate']*case['gas'] + P['onsite_vom_mwh']) return dict(engines=eng, fuel_vom=fuel, firm_gas=firm, battery=bess+losses, **b) def run(t=TARIFF, verbose=True): out={} for name, case in CASES.items(): A=config_A(case,t); C=config_C(case,t); Cc=config_C(case,t,credit_mw=100.0); E=config_E(case,t) tot=lambda d: sum(d.values()) r=dict(A=tot(A)/MWH, C=tot(C)/MWH, C_credit=tot(Cc)/MWH, E=tot(E)/MWH, diff=(tot(C)-tot(A))/MWH, diff_credit=(tot(Cc)-tot(A))/MWH, parts=dict(A=A,C=C,E=E)) out[name]=r if verbose: print(f"{name:24s} A ${r['A']:.1f} C ${r['C']:.1f} (credit ${r['C_credit']:.1f}) E ${r['E']:.1f} diff ${r['diff']:.1f} (credit ${r['diff_credit']:.1f})") return out if __name__=='__main__': if len(sys.argv)>1: j=json.load(open(sys.argv[1])) TARIFF.update({k:v for k,v in j.items() if k in TARIFF}) out=run() t=tariff_for(CASES['Base']) b=grid_bill(t) print(f"Base-case grid bill ${b/1e6:.1f}M = ${b/MWH:.1f}/MWh = {b/MWH/10:.2f} c/kWh") for name,r in out.items(): pa=r['parts']['A']; pc=r['parts']['C'] print(name, 'A parts $/MWh', {k:round(v/MWH,1) for k,v in pa.items()}) print(name, 'C parts $/MWh', {k:round(v/MWH,1) for k,v in pc.items()}) print(name, 'E parts $/MWh', {k:round(v/MWH,1) for k,v in r['parts']['E'].items()})