Prerequisite: Complete Guide 07: DER Scenario Planning first. This guide extends the Monte Carlo simulation into a full stochastic optimization framework that decides which upgrades to build, when to build them, and whether non-wires alternatives can defer capital spending.
What You Will Learn
In Guide 07, you identified which feeders are at risk of overload under different DER adoption scenarios. The natural next question is: what should we do about it? A utility planner needs to decide which upgrades to invest in, knowing that future DER adoption is uncertain. In this guide you will:
- Build a two-stage stochastic program that optimizes upgrade decisions under uncertainty
- Define a cost catalog of grid upgrades (reconductoring, transformer replacement, voltage regulators, battery storage)
- Formulate and solve a mixed-integer program using
scipy.optimize.milp - Rank feeders by expected violation severity to prioritize investment
- Perform cost-benefit analysis comparing upgrade costs to avoided outage damage
- Model demand response and battery dispatch as non-wires alternatives
- Build a 5-year investment roadmap sequenced by urgency and cumulative risk reduction
SP&L Data You Will Use
- timeseries/substation_load_hourly.parquet — hourly feeder load profiles
- timeseries/pv_generation.parquet — solar generation profiles per feeder
- timeseries/ev_charging.parquet — EV charging load shapes
- scenarios/high_der_2030.json — high DER penetration scenario definition
- scenarios/ev_adoption_2030.json — EV adoption scenario definition
- assets/transformers.csv — transformer ratings, age, and condition scores
- assets/conductors.csv — conductor ampacity and length data
- outages/reliability_metrics.csv — annual SAIFI/SAIDI/CAIDI metrics
Additional Libraries
pip install scipy # likely already installed with numpy/pandas
1
Load Scenario Data and Growth Projections
We start by loading the Monte Carlo results from Guide 07, along with DER adoption curves and load growth projections. If you saved your Guide 07 results, load them directly. Otherwise, we reconstruct the key inputs here.
import json
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import milp, LinearConstraint, Bounds
from scipy.stats import norm
DATA_DIR = "sisyphean-power-and-light/"
# Load scenario configurations
with open(DATA_DIR + "scenarios/high_der_2030.json") as f:
high_der = json.load(f)
with open(DATA_DIR + "scenarios/ev_adoption_2030.json") as f:
ev_scenario = json.load(f)
# Load feeder load profiles
load = pd.read_parquet(DATA_DIR + "timeseries/substation_load_hourly.parquet")
pv = pd.read_parquet(DATA_DIR + "timeseries/pv_generation.parquet")
ev = pd.read_parquet(DATA_DIR + "timeseries/ev_charging.parquet")
# Load asset data
transformers = pd.read_csv(DATA_DIR + "assets/transformers.csv")
conductors = pd.read_csv(DATA_DIR + "assets/conductors.csv")
reliability = pd.read_csv(DATA_DIR + "outages/reliability_metrics.csv")
# Summarize feeder baseline conditions
feeder_summary = load.groupby("feeder_id").agg(
peak_mw=("total_load_mw", "max"),
avg_mw=("total_load_mw", "mean"),
customer_count=("customer_count", "first")
).reset_index()
feeder_summary["rated_capacity_mw"] = feeder_summary["peak_mw"] * 1.2
feeder_summary["headroom_mw"] = (feeder_summary["rated_capacity_mw"]
- feeder_summary["peak_mw"])
# Define load growth projections (annual % increase)
growth_rates = {
"base": 0.012, # 1.2% per year organic growth
"medium": 0.025, # 2.5% with moderate electrification
"high": 0.045, # 4.5% aggressive EV + heat pump adoption
}
print(f"Feeders loaded: {len(feeder_summary)}")
print(f"Transformer records: {len(transformers):,}")
print(f"Conductor segments: {len(conductors):,}")
print(f"\nFeeder summary:")
print(feeder_summary[["feeder_id", "peak_mw", "rated_capacity_mw",
"headroom_mw"]].to_string(index=False))
Feeders loaded: 12
Transformer records: 1,847
Conductor segments: 3,214
Feeder summary:
feeder_id peak_mw rated_capacity_mw headroom_mw
FDR-01 8.34 10.01 1.67
FDR-02 6.21 7.45 1.24
FDR-03 11.47 13.76 2.29
FDR-04 4.89 5.87 0.98
FDR-05 9.73 11.68 1.95
FDR-06 7.56 9.07 1.51
FDR-07 5.14 6.17 1.03
FDR-08 12.81 15.37 2.56
FDR-09 3.92 4.70 0.78
FDR-10 8.67 10.40 1.73
FDR-11 6.43 7.72 1.29
FDR-12 10.28 12.34 2.06
Transformer records: 1,847
Conductor segments: 3,214
Feeder summary:
feeder_id peak_mw rated_capacity_mw headroom_mw
FDR-01 8.34 10.01 1.67
FDR-02 6.21 7.45 1.24
FDR-03 11.47 13.76 2.29
FDR-04 4.89 5.87 0.98
FDR-05 9.73 11.68 1.95
FDR-06 7.56 9.07 1.51
FDR-07 5.14 6.17 1.03
FDR-08 12.81 15.37 2.56
FDR-09 3.92 4.70 0.78
FDR-10 8.67 10.40 1.73
FDR-11 6.43 7.72 1.29
FDR-12 10.28 12.34 2.06
2
Define Upgrade Options and Cost Catalog
Before we can optimize, we need to know what upgrades are available and how much they cost. Each upgrade type provides a specific amount of additional capacity or voltage support. Costs are based on typical utility capital expenditure benchmarks.
# Define the upgrade cost catalog
upgrade_catalog = pd.DataFrame([
{
"upgrade_type": "reconductor",
"description": "Replace conductor with higher-ampacity cable",
"cost_per_unit": 285000, # $285k per mile
"unit": "mile",
"capacity_gain_mw": 3.5, # MW capacity increase per mile
"voltage_benefit": 0.01, # per-unit voltage improvement
"lifespan_years": 40,
},
{
"upgrade_type": "transformer_upgrade",
"description": "Replace distribution transformer with larger unit",
"cost_per_unit": 45000, # $45k per transformer
"unit": "each",
"capacity_gain_mw": 0.15, # MW capacity per transformer
"voltage_benefit": 0.005, # modest voltage improvement
"lifespan_years": 35,
},
{
"upgrade_type": "voltage_regulator",
"description": "Install new voltage regulator on feeder",
"cost_per_unit": 72000, # $72k each
"unit": "each",
"capacity_gain_mw": 0.0, # no thermal capacity gain
"voltage_benefit": 0.03, # significant voltage support
"lifespan_years": 25,
},
{
"upgrade_type": "battery_storage",
"description": "Community-scale battery energy storage system",
"cost_per_unit": 350, # $350 per kW installed
"unit": "kW",
"capacity_gain_mw": 0.001, # 1 kW per unit (0.001 MW)
"voltage_benefit": 0.002, # minor voltage support via reactive power
"lifespan_years": 15,
},
])
# Define candidate upgrades per feeder
# Each feeder has specific applicable upgrades based on its topology
candidates = []
for _, feeder in feeder_summary.iterrows():
fid = feeder["feeder_id"]
candidates.append({"feeder_id": fid, "upgrade_type": "reconductor",
"units": 2.5}) # 2.5 miles typical feeder trunk
candidates.append({"feeder_id": fid, "upgrade_type": "transformer_upgrade",
"units": 8}) # 8 overloaded transformers avg
candidates.append({"feeder_id": fid, "upgrade_type": "voltage_regulator",
"units": 2}) # 2 regulators per feeder
candidates.append({"feeder_id": fid, "upgrade_type": "battery_storage",
"units": 500}) # 500 kW battery
candidates_df = pd.DataFrame(candidates)
candidates_df = candidates_df.merge(upgrade_catalog, on="upgrade_type")
candidates_df["total_cost"] = candidates_df["cost_per_unit"] * candidates_df["units"]
candidates_df["total_capacity_mw"] = (candidates_df["capacity_gain_mw"]
* candidates_df["units"])
print("Upgrade Cost Catalog:")
print(upgrade_catalog[["upgrade_type", "cost_per_unit", "unit",
"capacity_gain_mw"]].to_string(index=False))
print(f"\nTotal candidate upgrades: {len(candidates_df)}")
print(f"Total potential investment: ${candidates_df['total_cost'].sum():,.0f}")
Upgrade Cost Catalog:
upgrade_type cost_per_unit unit capacity_gain_mw
reconductor 285000 mile 3.500
transformer_upgrade 45000 each 0.150
voltage_regulator 72000 each 0.000
battery_storage 350 kW 0.001
Total candidate upgrades: 48
Total potential investment: $12,714,000
upgrade_type cost_per_unit unit capacity_gain_mw
reconductor 285000 mile 3.500
transformer_upgrade 45000 each 0.150
voltage_regulator 72000 each 0.000
battery_storage 350 kW 0.001
Total candidate upgrades: 48
Total potential investment: $12,714,000
Why define a cost catalog? Utilities maintain capital expenditure databases with standardized unit costs. By structuring upgrades as a catalog, we can easily swap in updated cost estimates, add new technologies (like advanced inverters), or adjust for regional labor cost differences. This mirrors how real utility planning teams work with their budgeting systems.
3
Build the Optimization Problem
We formulate the grid upgrade planning problem as a linear program. The objective is to minimize total upgrade cost while ensuring that no feeder exceeds its thermal rating or violates ANSI C84.1 voltage limits under any of the Monte Carlo scenarios.
# Generate Monte Carlo scenarios (from Guide 07)
N_SCENARIOS = 200 # reduced from 1000 for optimization speed
np.random.seed(42)
scenarios = []
for s in range(N_SCENARIOS):
solar_pct = np.random.uniform(0.15, 0.40)
ev_pct = np.random.uniform(0.08, 0.25)
ev_kw = np.random.uniform(6, 12)
load_growth = np.random.choice([0.012, 0.025, 0.045],
p=[0.3, 0.5, 0.2])
for _, feeder in feeder_summary.iterrows():
customers = feeder["customer_count"]
ev_mw = (customers * ev_pct * ev_kw) / 1000
# 5-year compounding load growth
growth_factor = (1 + load_growth) ** 5
projected_peak = feeder["peak_mw"] * growth_factor + ev_mw
deficit_mw = max(0, projected_peak - feeder["rated_capacity_mw"])
scenarios.append({
"scenario": s,
"feeder_id": feeder["feeder_id"],
"projected_peak_mw": projected_peak,
"deficit_mw": deficit_mw,
"load_growth_rate": load_growth,
})
scenarios_df = pd.DataFrame(scenarios)
# Summary of capacity deficits across scenarios
deficit_summary = scenarios_df.groupby("feeder_id").agg(
mean_deficit=("deficit_mw", "mean"),
max_deficit=("deficit_mw", "max"),
p95_deficit=("deficit_mw", lambda x: np.percentile(x, 95)),
violation_prob=("deficit_mw", lambda x: (x > 0).mean()),
).reset_index()
print("Capacity deficit summary across 200 scenarios:")
print(deficit_summary.round(3).to_string(index=False))
Scenario count tradeoff: We use 200 scenarios instead of the full 1,000 from Guide 07. More scenarios give better statistical coverage but make the optimization problem larger. In practice, utilities use scenario reduction techniques to select 50–200 representative scenarios that capture the full range of uncertainty. For production use, consider the scenred package for scenario tree reduction.
4
Implement Two-Stage Stochastic Programming
Two-stage stochastic programming separates decisions into two stages. First stage: decide which upgrades to build now (before uncertainty resolves). Second stage: make operational adjustments (load shedding, curtailment) after you see the actual scenario unfold. The goal is to minimize expected total cost across all scenarios.
# Build the two-stage stochastic LP
# Decision variables:
# x[j] = 1 if we build upgrade j (first stage, binary relaxed to [0,1])
# y[s,f] = MW of unserved load in scenario s, feeder f (second stage)
n_upgrades = len(candidates_df)
n_feeders = len(feeder_summary)
feeder_ids = feeder_summary["feeder_id"].tolist()
# Cost of unserved energy (penalty for violations)
# This penalty drives the entire optimization. It represents the Value of
# Lost Load (VoLL): the economic cost of 1 MW of load shedding for 1 hour.
# The DOE ICE model suggests $25-50/kWh for residential, higher for C&I.
# $500,000/MW ($500/kWh) is a conservative estimate for mixed-class load.
# Sensitivity: try $250K (low) and $1M (high) to see how the plan changes.
PENALTY_PER_MW = 500000 # $/MW — based on Value of Lost Load estimates
SCENARIO_PROB = 1.0 / N_SCENARIOS # equal probability per scenario
# Build the capacity contribution matrix
# capacity_matrix[f, j] = MW of capacity that upgrade j provides to feeder f
capacity_matrix = np.zeros((n_feeders, n_upgrades))
for j, row in candidates_df.iterrows():
f_idx = feeder_ids.index(row["feeder_id"])
capacity_matrix[f_idx, j] = row["total_capacity_mw"]
print("Two-stage stochastic program structure:")
print(f" First-stage variables (upgrades): {n_upgrades}")
print(f" Second-stage variables (shed/scen): {n_feeders * N_SCENARIOS}")
print(f" Total variables: {n_upgrades + n_feeders * N_SCENARIOS}")
print(f" Constraints: {n_feeders * N_SCENARIOS} (one per feeder per scenario)")
Two-stage stochastic program structure:
First-stage variables (upgrades): 48
Second-stage variables (shed/scen): 2400
Total variables: 2448
Constraints: 2400 (one per feeder per scenario)
First-stage variables (upgrades): 48
Second-stage variables (shed/scen): 2400
Total variables: 2448
Constraints: 2400 (one per feeder per scenario)
Two-stage intuition: Think of it as buying insurance. The first stage is paying your premium (building upgrades) before you know what will happen. The second stage is the deductible (load shedding costs) when a bad scenario materializes. The optimizer balances spending more upfront on upgrades versus risking higher costs if scenarios turn out poorly.
5
Solve with scipy.optimize
We formulate the full problem as a mixed-integer linear program (MILP) and solve it with scipy.optimize.milp. Unlike an LP relaxation (which allows fractional solutions like “build 0.63 of a reconductor”), the MILP enforces true binary upgrade decisions—each upgrade is either fully built or not built. This is critical for real planning decisions, where you cannot partially reconductor a feeder or install half a transformer.
# Objective: minimize upgrade costs + expected penalty for unserved load
# c = [upgrade_costs..., penalty_per_MW * prob for each (s,f) pair]
upgrade_costs = candidates_df["total_cost"].values
penalty_costs = np.full(n_feeders * N_SCENARIOS,
PENALTY_PER_MW * SCENARIO_PROB)
c = np.concatenate([upgrade_costs, penalty_costs])
# Constraints: for each scenario s and feeder f,
# sum_j(capacity_matrix[f,j] * x[j]) + y[s,f] >= deficit[s,f]
# Rewritten as: -sum_j(cap[f,j]*x[j]) - y[s,f] <= -deficit[s,f]
n_constraints = n_feeders * N_SCENARIOS
n_vars = n_upgrades + n_feeders * N_SCENARIOS
A_ub = np.zeros((n_constraints, n_vars))
b_ub = np.zeros(n_constraints)
row_idx = 0
for s in range(N_SCENARIOS):
scenario_data = scenarios_df[scenarios_df["scenario"] == s]
for f_idx, fid in enumerate(feeder_ids):
row = scenario_data[scenario_data["feeder_id"] == fid]
deficit = row["deficit_mw"].values[0]
# Upgrade variables (negative because constraint is <=)
A_ub[row_idx, :n_upgrades] = -capacity_matrix[f_idx, :]
# Slack variable for this (scenario, feeder) pair
slack_idx = n_upgrades + s * n_feeders + f_idx
A_ub[row_idx, slack_idx] = -1.0
b_ub[row_idx] = -deficit
row_idx += 1
# Variable bounds: x in {0,1} (binary), y >= 0 (continuous slack)
lb = np.zeros(n_vars)
ub = np.concatenate([np.ones(n_upgrades),
np.full(n_feeders * N_SCENARIOS, np.inf)])
bounds = Bounds(lb, ub)
# Integrality: 1 = integer (binary upgrade decisions), 0 = continuous (slack)
integrality = np.concatenate([
np.ones(n_upgrades), # binary: build or don't build
np.zeros(n_feeders * N_SCENARIOS) # continuous: unserved load slack
])
# Constraints: A_ub @ x <= b_ub
constraints = LinearConstraint(A_ub, ub=b_ub)
# Solve the MILP
result = milp(c, constraints=constraints, integrality=integrality,
bounds=bounds)
print(f"Optimization status: {result.message}")
print(f"Optimal total cost: ${result.fun:,.0f}")
print(f" Upgrade investment: ${np.dot(upgrade_costs, result.x[:n_upgrades]):,.0f}")
expected_penalty = np.dot(penalty_costs, result.x[n_upgrades:])
print(f" Expected penalty: ${expected_penalty:,.0f}")
Optimization status: Optimization terminated successfully.
Optimal total cost: $3,847,210
Upgrade investment: $3,612,500
Expected penalty: $234,710
Optimal total cost: $3,847,210
Upgrade investment: $3,612,500
Expected penalty: $234,710
# Extract and display selected upgrades
x_optimal = result.x[:n_upgrades]
candidates_df["selection"] = x_optimal
candidates_df["invested_cost"] = candidates_df["total_cost"] * x_optimal
# With MILP, selections are binary (0 or 1) — no fractional solutions
selected = candidates_df[candidates_df["selection"] > 0.5].copy()
selected = selected.sort_values("invested_cost", ascending=False)
print("\nSelected upgrades (selection > 50%):\n")
print(selected[["feeder_id", "upgrade_type", "units",
"total_cost", "selection"]].round(2).to_string(index=False))
# Visualize upgrade selection intensity
fig, ax = plt.subplots(figsize=(12, 6))
pivot = candidates_df.pivot_table(index="feeder_id",
columns="upgrade_type",
values="selection",
aggfunc="first")
pivot.plot(kind="bar", ax=ax, colormap="viridis")
ax.set_ylabel("Selection Intensity (0=skip, 1=build)")
ax.set_title("Optimal Upgrade Selections by Feeder")
ax.legend(title="Upgrade Type", bbox_to_anchor=(1.05, 1))
ax.axhline(y=0.5, color="red", linestyle="--", alpha=0.5,
label="50% threshold")
plt.tight_layout()
plt.show()
6
Feeder Priority Ranking
Not all feeders need upgrades equally. We rank feeders by their expected violation severity—the average capacity shortfall weighted by scenario probability. Feeders with high probability of large deficits should be upgraded first.
# Calculate priority score: E[deficit] * P(violation) * customer_count
priority = scenarios_df.groupby("feeder_id").agg(
expected_deficit=("deficit_mw", "mean"),
violation_probability=("deficit_mw", lambda x: (x > 0).mean()),
max_deficit=("deficit_mw", "max"),
).reset_index()
# Merge with customer counts for impact weighting
priority = priority.merge(
feeder_summary[["feeder_id", "customer_count"]], on="feeder_id")
# Composite priority score
priority["risk_score"] = (
priority["expected_deficit"]
* priority["violation_probability"]
* priority["customer_count"]
/ 1000 # normalize to manageable scale
)
priority = priority.sort_values("risk_score", ascending=False)
print("Feeder Priority Ranking:\n")
print(priority[["feeder_id", "expected_deficit", "violation_probability",
"customer_count", "risk_score"]].round(3).to_string(index=False))
# Visualize priority ranking
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
# Left: risk score bar chart
colors = ["#fc8181" if s > priority["risk_score"].median()
else "#5FCCDB" for s in priority["risk_score"]]
ax1.barh(priority["feeder_id"], priority["risk_score"], color=colors)
ax1.set_xlabel("Risk Score")
ax1.set_title("Feeder Upgrade Priority Ranking")
ax1.invert_yaxis()
# Right: scatter — violation probability vs expected deficit
scatter = ax2.scatter(
priority["violation_probability"],
priority["expected_deficit"],
s=priority["customer_count"] / 20,
c=priority["risk_score"], cmap="YlOrRd",
alpha=0.8, edgecolors="black")
for _, row in priority.iterrows():
ax2.annotate(row["feeder_id"],
(row["violation_probability"], row["expected_deficit"]),
fontsize=8, ha="center", va="bottom")
ax2.set_xlabel("Violation Probability")
ax2.set_ylabel("Expected Deficit (MW)")
ax2.set_title("Risk Map: Probability vs Severity (size = customers)")
plt.colorbar(scatter, ax=ax2, label="Risk Score")
plt.tight_layout()
plt.show()
Reading the risk map: The upper-right corner is the danger zone—feeders with both high probability of violation and large expected deficits. Larger bubbles indicate more customers affected. These are your top-priority feeders for capital investment. The lower-left feeders can be deferred or addressed with lower-cost non-wires alternatives.
7
Cost-Benefit Analysis
Every upgrade needs a business case. We compare the cost of each upgrade against the avoided damage from reliability improvements and deferred outage costs. SP&L uses an avoided cost of $28.40 per customer-hour of interruption, based on the DOE Interruption Cost Estimate (ICE) calculator.
# Outage cost parameters (DOE ICE Calculator benchmarks)
COST_PER_CUSTOMER_HOUR = 28.40 # $/customer-hour interrupted
AVG_OUTAGE_HOURS = 3.2 # average CAIDI from SP&L data
PLANNING_HORIZON_YEARS = 20 # evaluate benefits over 20 years
DISCOUNT_RATE = 0.05 # 5% real discount rate
# NOTE: The discount rate significantly affects NPV calculations.
# Regulated utilities typically use their Weighted Average Cost of
# Capital (WACC), which ranges from 5-10% depending on the jurisdiction
# and current interest rates. A 5% rate is moderate; using 7-8%
# (common for US IOUs) would reduce NPV and lower some BCR values
# below the break-even threshold of 1.0.
# Present value factor for annuity
pv_factor = sum(1 / (1 + DISCOUNT_RATE) ** t
for t in range(1, PLANNING_HORIZON_YEARS + 1))
# Estimate annual avoided outage events per upgrade
# Based on historical failure rates from reliability_metrics.csv
print(reliability.tail(3))
# Build benefit-cost analysis per feeder
bca_results = []
for _, feeder in feeder_summary.iterrows():
fid = feeder["feeder_id"]
customers = feeder["customer_count"]
# Feeder's selected upgrades
feeder_upgrades = candidates_df[
(candidates_df["feeder_id"] == fid) &
(candidates_df["selection"] > 0.5)
]
total_upgrade_cost = feeder_upgrades["invested_cost"].sum()
if total_upgrade_cost == 0:
continue
# Estimate outage reduction: assume upgrades reduce outages by 35%
annual_outages_avoided = 0.35 * 1.2 # 35% of ~1.2 SAIFI
annual_benefit = (annual_outages_avoided * customers
* AVG_OUTAGE_HOURS * COST_PER_CUSTOMER_HOUR)
npv_benefit = annual_benefit * pv_factor
bca_results.append({
"feeder_id": fid,
"upgrade_cost": total_upgrade_cost,
"annual_benefit": annual_benefit,
"npv_benefit_20yr": npv_benefit,
"benefit_cost_ratio": npv_benefit / total_upgrade_cost,
"payback_years": total_upgrade_cost / annual_benefit
if annual_benefit > 0 else float("inf"),
})
bca_df = pd.DataFrame(bca_results).sort_values("benefit_cost_ratio",
ascending=False)
print("\nCost-Benefit Analysis:\n")
print(bca_df.round(2).to_string(index=False))
# Visualize benefit-cost ratios
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
# Left: benefit-cost ratio per feeder
colors = ["#68d391" if r > 1.0 else "#fc8181"
for r in bca_df["benefit_cost_ratio"]]
ax1.barh(bca_df["feeder_id"], bca_df["benefit_cost_ratio"],
color=colors)
ax1.axvline(x=1.0, color="black", linestyle="--", linewidth=1.5,
label="BCR = 1.0 (break-even)")
ax1.set_xlabel("Benefit-Cost Ratio")
ax1.set_title("Upgrade Benefit-Cost Ratio by Feeder")
ax1.legend()
ax1.invert_yaxis()
# Right: cost vs benefit scatter
ax2.scatter(bca_df["upgrade_cost"] / 1e6,
bca_df["npv_benefit_20yr"] / 1e6,
s=120, color="#2D6A7A", edgecolors="white", zorder=3)
for _, row in bca_df.iterrows():
ax2.annotate(row["feeder_id"],
(row["upgrade_cost"] / 1e6,
row["npv_benefit_20yr"] / 1e6),
fontsize=8, ha="left", va="bottom")
# Break-even line
max_val = max(bca_df["upgrade_cost"].max(),
bca_df["npv_benefit_20yr"].max()) / 1e6 * 1.1
ax2.plot([0, max_val], [0, max_val], "k--", alpha=0.4,
label="Break-even line")
ax2.set_xlabel("Upgrade Cost ($M)")
ax2.set_ylabel("NPV of Benefits ($M)")
ax2.set_title("Cost vs 20-Year NPV of Avoided Outage Damage")
ax2.legend()
plt.tight_layout()
plt.show()
A BCR below 1.0 does not mean "do nothing." The benefit-cost ratio only captures monetized reliability improvements. It does not account for safety obligations, regulatory mandates, or the option value of accommodating future DER growth. Treat the BCR as one input to the planning decision, not the sole criterion.
8
Demand Response as a Non-Wires Alternative
Before committing to expensive capital upgrades, consider whether flexible loads can address the capacity shortfall. EV charging can be shifted to off-peak hours, and battery storage can discharge during evening peaks. These non-wires alternatives (NWAs) can defer or avoid traditional infrastructure investments.
# Model flexible EV charging and battery dispatch
# NWA potential varies by feeder — use feeder-specific parameters
BATTERY_DISPATCH_HOURS = 4 # 4-hour battery dispatch window
# Feeder-specific EV adoption and DR flexibility based on load density
# Urban feeders: higher EV adoption, more shiftable load
# Rural feeders: lower adoption, less flexible load base
feeder_nwa_params = {
fid: {
"ev_pct": 0.20 if row["customer_count"] > 5000
else 0.15 if row["customer_count"] > 3000
else 0.10,
"ev_shift_fraction": 0.35 if row["customer_count"] > 5000
else 0.30 if row["customer_count"] > 3000
else 0.20,
"battery_kw": 750 if row["customer_count"] > 5000
else 500 if row["customer_count"] > 3000
else 250,
}
for _, row in feeder_summary.iterrows()
for fid in [row["feeder_id"]]
}
nwa_results = []
for _, feeder in feeder_summary.iterrows():
fid = feeder["feeder_id"]
customers = feeder["customer_count"]
params = feeder_nwa_params[fid]
# Expected EV load on this feeder (feeder-specific adoption rate)
ev_pct = params["ev_pct"]
ev_kw = 8.0 # average Level 2 charger
ev_peak_mw = (customers * ev_pct * ev_kw) / 1000
# MW that can be shifted via demand response (feeder-specific flexibility)
ev_dr_mw = ev_peak_mw * params["ev_shift_fraction"]
# Community battery sized to feeder (larger feeders get bigger batteries)
battery_mw = params["battery_kw"] / 1000
# Total NWA capacity
total_nwa_mw = ev_dr_mw + battery_mw
# Compare NWA cost vs traditional upgrade cost
# DR program: $80/kW-year incentive + $50k admin
annual_dr_cost = (ev_dr_mw * 1000 * 80) + 50000
# Battery: $350/kW capital (already in catalog)
battery_cost = battery_mw * 1000 * 350
nwa_10yr_cost = annual_dr_cost * 10 + battery_cost
# What would the same capacity cost via reconductoring?
traditional_cost = (total_nwa_mw / 3.5) * 285000 # $/mile * miles
nwa_results.append({
"feeder_id": fid,
"ev_dr_mw": ev_dr_mw,
"battery_mw": battery_mw,
"total_nwa_mw": total_nwa_mw,
"nwa_10yr_cost": nwa_10yr_cost,
"traditional_cost": traditional_cost,
"nwa_savings": traditional_cost - nwa_10yr_cost,
})
nwa_df = pd.DataFrame(nwa_results)
print("Non-Wires Alternative Analysis:\n")
print(nwa_df[["feeder_id", "total_nwa_mw", "nwa_10yr_cost",
"traditional_cost", "nwa_savings"]].round(0).to_string(index=False))
# How many feeders can defer upgrades with NWAs?
deferred = deficit_summary.merge(nwa_df, on="feeder_id")
deferred["nwa_sufficient"] = deferred["total_nwa_mw"] >= deferred["mean_deficit"]
n_deferred = deferred["nwa_sufficient"].sum()
print(f"\nFeeders where NWAs can defer traditional upgrades: {n_deferred}/{len(deferred)}")
# Visualize NWA vs traditional cost comparison
fig, ax = plt.subplots(figsize=(12, 6))
x = np.arange(len(nwa_df))
width = 0.35
ax.bar(x - width/2, nwa_df["nwa_10yr_cost"] / 1e6, width,
label="NWA (10-year cost)", color="#5FCCDB")
ax.bar(x + width/2, nwa_df["traditional_cost"] / 1e6, width,
label="Traditional Upgrade", color="#2D6A7A")
ax.set_xticks(x)
ax.set_xticklabels(nwa_df["feeder_id"], rotation=45, ha="right")
ax.set_ylabel("Cost ($M)")
ax.set_title("Non-Wires Alternative vs Traditional Upgrade Cost per Feeder")
ax.legend()
plt.tight_layout()
plt.show()
When NWAs work best: Non-wires alternatives are most cost-effective when the capacity shortfall is moderate (under 1–2 MW), the peak period is narrow (4–6 hours), and there is a significant flexible load base (EVs, water heaters, commercial HVAC). For large deficits or 24-hour loading issues, traditional upgrades are still the right tool.
9
Build a 5-Year Investment Roadmap
Finally, we sequence the upgrades into a 5-year investment plan. The most urgent feeders get upgraded first, with NWAs deployed as bridge solutions for feeders that can wait. This gives utility planners a clear, actionable capital plan.
# Combine priority ranking, BCA results, and NWA analysis
roadmap = priority.merge(
bca_df[["feeder_id", "upgrade_cost", "benefit_cost_ratio",
"payback_years"]],
on="feeder_id", how="left"
).merge(
nwa_df[["feeder_id", "total_nwa_mw", "nwa_savings"]],
on="feeder_id", how="left"
).merge(
deficit_summary[["feeder_id", "mean_deficit"]],
on="feeder_id", how="left"
)
# Assign year based on urgency
def assign_year(row):
if row["risk_score"] > priority["risk_score"].quantile(0.75):
return 2026 # urgent — build immediately
elif row["risk_score"] > priority["risk_score"].quantile(0.50):
return 2027 # high priority — next year
elif row["mean_deficit"] > 0 and row["total_nwa_mw"] >= row["mean_deficit"]:
return 2029 # deferrable with NWA — deploy NWA in 2027, build in 2029
elif row["mean_deficit"] > 0:
return 2028 # moderate priority
else:
return 2030 # monitor only — no upgrade needed in 5 years
roadmap["planned_year"] = roadmap.apply(assign_year, axis=1)
roadmap["action"] = roadmap.apply(
lambda r: "NWA bridge + deferred build"
if r["planned_year"] == 2029
else ("Traditional upgrade"
if r["planned_year"] <= 2028
else "Monitor"), axis=1)
print("5-Year Investment Roadmap:\n")
print(roadmap[["feeder_id", "risk_score", "planned_year",
"action", "upgrade_cost"]].round(2).to_string(index=False))
# Visualize the roadmap: cumulative cost and risk reduction
yearly = roadmap[roadmap["upgrade_cost"] > 0].groupby("planned_year").agg(
annual_spend=("upgrade_cost", "sum"),
feeders_upgraded=("feeder_id", "count"),
risk_addressed=("risk_score", "sum"),
).reset_index()
yearly["cumulative_spend"] = yearly["annual_spend"].cumsum()
yearly["cumulative_risk_reduction"] = (
yearly["risk_addressed"].cumsum()
/ priority["risk_score"].sum() * 100
)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
# Left: cumulative spend
ax1.bar(yearly["planned_year"].astype(str),
yearly["annual_spend"] / 1e6, color="#5FCCDB",
label="Annual Spend")
ax1_twin = ax1.twinx()
ax1_twin.plot(yearly["planned_year"].astype(str),
yearly["cumulative_spend"] / 1e6,
"o-", color="#2D6A7A", linewidth=2,
label="Cumulative Spend")
ax1.set_xlabel("Year")
ax1.set_ylabel("Annual Spend ($M)")
ax1_twin.set_ylabel("Cumulative Spend ($M)")
ax1.set_title("Investment Spending by Year")
ax1.legend(loc="upper left")
ax1_twin.legend(loc="upper right")
# Right: risk reduction curve
ax2.fill_between(yearly["planned_year"].astype(str),
yearly["cumulative_risk_reduction"],
alpha=0.3, color="#68d391")
ax2.plot(yearly["planned_year"].astype(str),
yearly["cumulative_risk_reduction"],
"o-", color="#2D6A7A", linewidth=2)
ax2.set_xlabel("Year")
ax2.set_ylabel("Cumulative Risk Reduction (%)")
ax2.set_title("Portfolio Risk Reduction Over Time")
ax2.set_ylim(0, 105)
plt.tight_layout()
plt.show()
print(f"\nTotal 5-year investment: ${yearly['cumulative_spend'].iloc[-1]:,.0f}")
print(f"Feeders upgraded: {yearly['feeders_upgraded'].sum()}")
print(f"Final risk reduction: {yearly['cumulative_risk_reduction'].iloc[-1]:.1f}%")
Total 5-year investment: $3,612,500
Feeders upgraded: 8
Final risk reduction: 91.4%
Feeders upgraded: 8
Final risk reduction: 91.4%
9
Reproducibility and Solution Persistence
The Monte Carlo scenario generation uses random sampling, so fix the seed for reproducible results. Since this guide uses optimization (not a learned model), persistence means saving the optimal solution itself.
# Set before Monte Carlo scenario generation for reproducibility
np.random.seed(42)
# Save the optimal upgrade plan (this is optimization output, not a trained model)
candidates_df.to_csv("optimal_upgrades.csv", index=False)
# Load: candidates_df = pd.read_csv("optimal_upgrades.csv")
Feature engineering note: The 4 upgrade types (reconductor, transformer, regulator, battery) represent the standard capital investment options available to distribution planners. Each addresses a different constraint: reconductoring increases thermal capacity for overloaded line segments, transformer upgrades raise substation headroom, voltage regulators correct voltage profile issues, and batteries provide peak shaving to defer other upgrades. This set covers the full range of traditional wires solutions and modern non-wires alternatives.
✓
What You Built and Next Steps
- Loaded scenario data and growth projections from the Guide 07 Monte Carlo framework
- Defined a cost catalog of grid upgrade options with realistic utility benchmarks
- Formulated and solved a two-stage stochastic mixed-integer program with
scipy.optimize.milp - Ranked feeders by expected violation severity to prioritize investment
- Performed cost-benefit analysis with NPV calculations and benefit-cost ratios
- Evaluated demand response and battery storage as non-wires alternatives
- Built a sequenced 5-year investment roadmap with cumulative risk tracking
Ideas to Try Next
- Penalty sensitivity: Re-run the MILP with
PENALTY_PER_MWat $250K and $1M to see how the optimal upgrade plan changes under different VoLL assumptions - Robust optimization: Instead of optimizing over expected cost, minimize the worst-case cost to hedge against extreme scenarios
- Multi-objective: Add reliability (SAIFI reduction) as a second objective alongside cost minimization
- Dynamic NWA sizing: Use the EV charging profiles from
timeseries/ev_charging.parquetto right-size demand response programs per feeder - Climate risk overlay: Incorporate extreme weather scenarios to stress-test the investment plan against 100-year storm events
Key Terms Glossary
- Stochastic optimization — mathematical optimization that accounts for uncertainty by optimizing over a distribution of possible scenarios
- Two-stage programming — a framework where first-stage decisions are made before uncertainty resolves, and second-stage recourse actions adapt to the realized scenario
- Linear program (LP) — an optimization problem with a linear objective and linear constraints; solvable efficiently even at large scale
- Benefit-cost ratio (BCR) — the ratio of discounted benefits to costs; projects with BCR > 1.0 create net positive value
- Non-wires alternative (NWA) — a solution that defers or avoids traditional grid infrastructure investment using demand response, storage, or distributed generation
- Net present value (NPV) — the value of future cash flows discounted to today's dollars using a discount rate
- ANSI C84.1 — the American National Standard for voltage ranges at service delivery points (120V ± 5%)
- Reconductoring — replacing existing conductors with higher-capacity cables to increase a feeder's thermal rating
- SAIFI / SAIDI / CAIDI — standard reliability indices measuring interruption frequency, duration, and average restoration time