Heat Smart Orkney: Quantifying Wind Energy Curtailment and the Commercial Case for Demand Response¶
Group L¶
Note:¶
This notebook is not the appendix section of the accompanying report. This notebook maps directly to the appendix sections of the accompanying report, which contains all narrative, methodology, assumptions, figures, tables, results and conclusions. Each cell corresponds to a named appendix section and contains the underlying code used to produce the outputs documented there. This notebook is provided solely for reproducibility.¶
Data Loading and Clean¶
AppendixDataLoadAndClean.1 — Imports and Constants¶
In [28]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline, interp1d
import matplotlib as mpl
# Turbine constants (Assumptions Register v5 §2)
RATED_CAPACITY_KW = 900.0
CUT_IN_MS = 2.5
STORM_TAPER_START = 28.0
CUT_OUT_MS = 34.0
HOURS_PER_STEP = 1 / 60
BIN_WIDTH = 0.5
MIN_OBS = 50
DATA_PATH = 'data/Turbine_telemetry.csv'
AppendixDataLoadAndClean.2 Load Raw Data¶
In [29]:
# AppendixDataLoadAndClean.2 — Load Raw Data
df = pd.read_csv(DATA_PATH, parse_dates=['Timestamp'])
df = df.sort_values('Timestamp').reset_index(drop=True)
print(f"Rows loaded: {len(df):,}")
print(f"Date range: {df['Timestamp'].min()} → {df['Timestamp'].max()}")
print(f"\nMissing values:\n{df.isnull().sum()}")
print(f"\nBasic statistics:")
print(df[['Power_kw', 'Setpoint_kw', 'Wind_ms']].describe().round(2))
Rows loaded: 1,069,636
Date range: 2015-05-28 00:00:23 → 2018-01-11 06:14:32
Missing values:
Timestamp 0
Power_kw 2702
Setpoint_kw 2715
Wind_ms 2705
dtype: int64
Basic statistics:
Power_kw Setpoint_kw Wind_ms
count 1066934.00 1066921.00 1066931.00
mean 367.19 810.66 9.96
std 339.96 261.51 5.45
min 0.00 0.00 0.00
25% 51.00 900.00 6.10
50% 251.00 900.00 9.00
75% 704.00 900.00 12.40
max 938.00 900.00 47.20
AppendixDataLoadAndClean.3 — Drop Missing Values¶
In [30]:
n_before = len(df)
df = df.dropna(subset=['Power_kw', 'Setpoint_kw', 'Wind_ms']).copy()
df = df.reset_index(drop=True)
print(f"Rows before: {n_before:,}")
print(f"Rows dropped: {n_before - len(df):,} ({(n_before - len(df))/n_before*100:.2f}%)")
print(f"Rows remaining: {len(df):,}")
Rows before: 1,069,636 Rows dropped: 5,417 (0.51%) Rows remaining: 1,064,219
AppendixDataLoadAndClean.4 — Validate Zero Power Readings¶
In [31]:
# Zeros are physically valid: turbine below cut-in or halted by operator.
# Do NOT drop zeros.
zero_power = (df['Power_kw'] == 0)
zero_wind = (df['Wind_ms'] == 0)
print(f"Zero power rows: {zero_power.sum():,} ({zero_power.mean()*100:.1f}%)")
print(f"Zero wind rows: {zero_wind.sum():,} ({zero_wind.mean()*100:.1f}%)")
print(f"Zero power + non-zero wind: {(zero_power & ~zero_wind).sum():,} — likely stopped/curtailed")
Zero power rows: 130,760 (12.3%) Zero wind rows: 1,142 (0.1%) Zero power + non-zero wind: 129,618 — likely stopped/curtailed
AppendixDataLoadAndClean.5 — Setpoint Distribution¶
In [32]:
print("Setpoint distribution (top 5 values):")
sp_counts = df['Setpoint_kw'].value_counts().head(5)
for val, count in sp_counts.items():
print(f" {val:>6.0f} kW: {count:>8,} rows ({count/len(df)*100:.1f}%)")
Setpoint distribution (top 5 values):
900 kW: 943,075 rows (88.6%)
0 kW: 89,806 rows (8.4%)
189 kW: 1,470 rows (0.1%)
500 kW: 938 rows (0.1%)
549 kW: 511 rows (0.0%)
AppendixDataLoadAndClean.6 — Assign Operational State¶
In [33]:
# spinning = setpoint 900 kW, normal operation
# curtailed = setpoint between 0 and 900 (exclusive), operator-constrained
# stopped = setpoint 0 kW, turbine halted
def assign_state(row):
if row['Setpoint_kw'] == RATED_CAPACITY_KW:
return 'spinning'
elif row['Setpoint_kw'] == 0:
return 'stopped'
else:
return 'curtailed'
df['state'] = df.apply(assign_state, axis=1)
print("State distribution:")
for state, count in df['state'].value_counts().items():
print(f" {state:>12}: {count:>8,} rows ({count/len(df)*100:.1f}%)")
State distribution:
spinning: 943,075 rows (88.6%)
stopped: 89,806 rows (8.4%)
curtailed: 31,338 rows (2.9%)
AppendixDataLoadAndClean.7 — Data Coverage Plot¶
In [34]:
fig, ax = plt.subplots(figsize=(14, 3))
daily = df.set_index('Timestamp').resample('D').size()
ax.bar(daily.index, daily.values, width=1, color='steelblue', alpha=0.7)
ax.axhline(y=1440, color='red', linestyle='--', alpha=0.5, label='Full day (1440 min)')
ax.set_ylabel('Records per day')
ax.set_title('Figure A3.1: Daily Data Coverage (May 2015 – Jan 2018)')
ax.legend()
plt.tight_layout()
plt.show()
AppendixDataLoadAndClean.8 — Confirm Output¶
In [35]:
print("=== AppendixDataLoadAndClean complete ===")
print(f" Rows: {len(df):,}")
print(f" Columns: {list(df.columns)}")
print(f" States: {df['state'].value_counts().to_dict()}")
=== AppendixDataLoadAndClean complete ===
Rows: 1,064,219
Columns: ['Timestamp', 'Power_kw', 'Setpoint_kw', 'Wind_ms', 'state']
States: {'spinning': 943075, 'stopped': 89806, 'curtailed': 31338}
Theme and Colors¶
In [6]:
# ── HSO Report Style ──────────────────────────────────────────
# C1 #0B1020 near black — primary curves, fitted lines
# C2 #FF5A6E red/coral — curtailed energy, alerts
# C3 #12C2A9 teal/cyan — actual generation
# C4 #7C5CFF purple — potential energy
# C5 #243B53 dark navy — scatter observations, secondary
# Background fills: #F4F7FB
# ─────────────────────────────────────────────────────────────
C1 = '#0B1020'
C2 = '#FF5A6E'
C3 = '#12C2A9'
C4 = '#7C5CFF'
C5 = '#243B53'
BG = '#F4F7FB'
HSO_PALETTE = [C1, C2, C3, C4, C5]
mpl.rcParams.update({
'figure.facecolor': 'white',
'figure.dpi': 150,
'axes.facecolor': 'white',
'axes.edgecolor': '#CCCCCC',
'axes.linewidth': 0.8,
'axes.spines.top': False,
'axes.spines.right': False,
'axes.grid': True,
'grid.color': '#EEEEEE',
'grid.linewidth': 0.6,
'grid.linestyle': '-',
'xtick.color': '#444444',
'ytick.color': '#444444',
'xtick.labelsize': 10,
'ytick.labelsize': 10,
'xtick.major.size': 3,
'ytick.major.size': 3,
'axes.labelsize': 11,
'axes.labelcolor': '#111111',
'axes.labelweight': 'normal',
'axes.titlesize': 12,
'axes.titleweight': 'bold',
'axes.titlecolor': '#111111',
'axes.titlepad': 10,
'legend.frameon': False,
'legend.fontsize': 9,
'legend.labelcolor': '#111111',
'font.family': 'sans-serif',
'font.sans-serif': ['Arial', 'Helvetica Neue', 'DejaVu Sans'],
'text.color': '#111111',
})
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[6], line 19 15 BG = '#F4F7FB' 17 HSO_PALETTE = [C1, C2, C3, C4, C5] ---> 19 mpl.rcParams.update({ 20 'figure.facecolor': 'white', 21 'figure.dpi': 150, 22 'axes.facecolor': 'white', 23 'axes.edgecolor': '#CCCCCC', 24 'axes.linewidth': 0.8, 25 'axes.spines.top': False, 26 'axes.spines.right': False, 27 'axes.grid': True, 28 'grid.color': '#EEEEEE', 29 'grid.linewidth': 0.6, 30 'grid.linestyle': '-', 31 'xtick.color': '#444444', 32 'ytick.color': '#444444', 33 'xtick.labelsize': 10, 34 'ytick.labelsize': 10, 35 'xtick.major.size': 3, 36 'ytick.major.size': 3, 37 'axes.labelsize': 11, 38 'axes.labelcolor': '#111111', 39 'axes.labelweight': 'normal', 40 'axes.titlesize': 12, 41 'axes.titleweight': 'bold', 42 'axes.titlecolor': '#111111', 43 'axes.titlepad': 10, 44 'legend.frameon': False, 45 'legend.fontsize': 9, 46 'legend.labelcolor': '#111111', 47 'font.family': 'sans-serif', 48 'font.sans-serif': ['Arial', 'Helvetica Neue', 'DejaVu Sans'], 49 'text.color': '#111111', 50 }) NameError: name 'mpl' is not defined
Appendix B: Power Curve Construction and Validation¶
Appendix B.1 — Raw Scatter¶
In [165]:
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
state_colours = {'spinning': C5, 'curtailed': C2, 'stopped': C4}
for state, colour in state_colours.items():
mask = df['state'] == state
axes[0].scatter(df.loc[mask, 'Wind_ms'], df.loc[mask, 'Power_kw'],
s=0.3, alpha=0.15, c=colour, label=state.capitalize())
axes[0].axhline(y=RATED_CAPACITY_KW, color=C2, linestyle='--',
linewidth=0.8, label='Rated 900 kW')
axes[0].set_xlabel('Wind speed (m/s)')
axes[0].set_ylabel('Power output (kW)')
axes[0].set_title('Figure B.1: Raw Scatter — All Observations by State')
axes[0].legend(markerscale=6, fontsize=9)
axes[1].scatter(df['Wind_ms'], df['Power_kw'],
s=0.3, alpha=0.05, c=C5)
axes[1].axhline(y=RATED_CAPACITY_KW, color=C2, linestyle='--',
linewidth=0.8, label='Rated 900 kW')
axes[1].set_xlabel('Wind speed (m/s)')
axes[1].set_ylabel('Power output (kW)')
axes[1].set_title('Figure B.2: Raw Scatter — Density Adjusted')
axes[1].legend(fontsize=9)
plt.tight_layout()
plt.savefig('fig_5_1_2_raw_scatter.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"Total observations: {len(df):,}")
print(f"State breakdown:")
for state, count in df['state'].value_counts().items():
print(f" {state:>12}: {count:>8,} ({count/len(df)*100:.1f}%)")
Total observations: 1,064,219
State breakdown:
spinning: 943,075 (88.6%)
stopped: 89,806 (8.4%)
curtailed: 31,338 (2.9%)
Appendix B.2 — Filtering Unconstrained Records¶
In [38]:
# Unconstrained training sample:
# 1. Setpoint_kw = 900 — no explicit setpoint reduction
# 2. Power_kw > 0 — turbine generating
# 3. Wind_ms > 0 — wind resource present
# 4. Power_kw <= Setpoint_kw — removes transient overshoots that would
# bias the upper envelope above the 900 kW rated limit
uc = df[
(df['Setpoint_kw'] == RATED_CAPACITY_KW) &
(df['Power_kw'] > 0) &
(df['Wind_ms'] > 0) &
(df['Power_kw'] <= df['Setpoint_kw'])
].copy()
print(f"Full dataset: {len(df):,} rows")
print(f"Unconstrained sample: {len(uc):,} rows ({len(uc)/len(df)*100:.1f}%)")
print(f"Excluded (curtailed/stopped): {len(df) - len(uc):,} rows")
fig, ax = plt.subplots(figsize=(12, 5))
excluded = df[~df.index.isin(uc.index)]
ax.scatter(excluded['Wind_ms'], excluded['Power_kw'],
s=0.3, alpha=0.1, c=C4, label='Excluded (curtailed/stopped)')
ax.scatter(uc['Wind_ms'], uc['Power_kw'],
s=0.3, alpha=0.1, c=C3, label='Unconstrained sample')
ax.axhline(y=RATED_CAPACITY_KW, color=C2, linestyle='--',
linewidth=0.8, label='Rated 900 kW')
ax.set_xlabel('Wind speed (m/s)')
ax.set_ylabel('Power output (kW)')
ax.set_title('Figure 5.3: Unconstrained Sample vs Excluded Observations')
ax.legend(markerscale=6, fontsize=9)
plt.tight_layout()
plt.savefig('fig_5_3_filtered_scatter.png', dpi=150, bbox_inches='tight')
plt.show()
Full dataset: 1,064,219 rows Unconstrained sample: 785,075 rows (73.8%) Excluded (curtailed/stopped): 279,144 rows
Appendix B.3 — Binned P50 / P90 / P95¶
In [166]:
### AppendixPowerCurve.3 — Binning and P50/P90/P95 Estimation
uc['wind_bin'] = (uc['Wind_ms'] / BIN_WIDTH).round() * BIN_WIDTH
power_curve_df = uc.groupby('wind_bin').agg(
p50_power=('Power_kw', 'median'),
p90_power=('Power_kw', lambda x: x.quantile(0.90)),
p95_power=('Power_kw', lambda x: x.quantile(0.95)),
std_power=('Power_kw', 'std'),
count=('Power_kw', 'count')
).reset_index()
# Remove sparse bins
power_curve_df = power_curve_df[power_curve_df['count'] >= MIN_OBS].copy()
# Cap all estimators at rated capacity
for col in ['p50_power', 'p90_power', 'p95_power']:
power_curve_df[f'{col}_capped'] = power_curve_df[col].clip(upper=RATED_CAPACITY_KW)
print(f"Wind speed range: {power_curve_df['wind_bin'].min():.1f} – "
f"{power_curve_df['wind_bin'].max():.1f} m/s")
print(f"Number of bins (>={MIN_OBS} obs): {len(power_curve_df)}")
print(f"\n{'Wind':>6} {'P50':>8} {'P90':>8} {'P95':>8} {'Count':>8}")
for _, row in power_curve_df[
power_curve_df['wind_bin'].isin([3, 5, 7, 9, 11, 13, 15, 18, 22, 26])
].iterrows():
print(f" {row['wind_bin']:>4.1f} {row['p50_power_capped']:>8.0f} "
f"{row['p90_power_capped']:>8.0f} {row['p95_power_capped']:>8.0f} "
f"{row['count']:>8.0f}")
# Truncate both scatter and bin lines at 28 m/s
# Storm taper is an imposed physical constraint, not empirically estimated
# Showing sparse noisy bins above 28 m/s would undermine confidence in the estimators
uc_plot = uc[uc['Wind_ms'] < STORM_TAPER_START]
plot_df = power_curve_df[power_curve_df['wind_bin'] < STORM_TAPER_START]
fig, ax = plt.subplots(figsize=(12, 6))
ax.scatter(uc_plot['Wind_ms'], uc_plot['Power_kw'],
s=0.3, alpha=0.05, c=C5, label='Unconstrained observations')
ax.plot(plot_df['wind_bin'], plot_df['p50_power_capped'],
'o-', color=C4, markersize=3, linewidth=1.2, alpha=0.7, label='P50 (median)')
ax.plot(plot_df['wind_bin'], plot_df['p90_power_capped'],
's-', color=C3, markersize=4, linewidth=2, label='P90 (base case)')
ax.plot(plot_df['wind_bin'], plot_df['p95_power_capped'],
'^-', color=C1, markersize=4, linewidth=2, alpha=0.8, label='P95 (sensitivity)')
ax.axhline(y=RATED_CAPACITY_KW, color=C2, linestyle='--',
linewidth=1, label='Rated 900 kW')
ax.axvspan(STORM_TAPER_START, CUT_OUT_MS, color=C2, alpha=0.07)
ax.text(31, 50, 'Storm taper\n28→34 m/s', ha='center', fontsize=9, color=C2)
ax.set_xlabel('Wind speed (m/s)')
ax.set_ylabel('Power output (kW)')
ax.set_title('Figure B.3: Binned P50 / P90 / P95 over Unconstrained Scatter')
ax.legend(markerscale=4, fontsize=9)
ax.set_xlim(0, 28.5)
ax.set_ylim(-20, 1000)
plt.tight_layout()
plt.savefig('fig_5_4_binned_estimators.png', dpi=150, bbox_inches='tight')
plt.show()
Wind speed range: 2.5 – 35.0 m/s Number of bins (>=50 obs): 66 Wind P50 P90 P95 Count 3.0 5 10 12 16131 5.0 56 77 83 40515 7.0 170 209 224 42562 9.0 366 438 459 37951 11.0 657 748 773 40721 13.0 895 900 900 8445 15.0 896 900 900 4264 18.0 896 900 900 3131 22.0 895 899 900 2023 26.0 892 899 900 1025
Appendix B.4 — Spline Fitting and Physical Constraint Enforcement¶
In [40]:
from scipy.interpolate import UnivariateSpline, interp1d
# Smoothing factor s = n * 50
# Lower s → oscillation in 6-12 m/s region
# Higher s → loses cubic shape across rising region
# Verified visually: current value tracks bin centres without artefacts
curve_data = power_curve_df[
['wind_bin', 'p50_power_capped', 'p90_power_capped', 'p95_power_capped']
].copy()
splines = {}
spline_fit = {}
for label, col in [
('p50', 'p50_power_capped'),
('p90', 'p90_power_capped'),
('p95', 'p95_power_capped')
]:
try:
s_val = len(curve_data) * 50
sp = UnivariateSpline(curve_data['wind_bin'],
curve_data[col], s=s_val, k=3)
splines[label] = sp
spline_fit[label] = True
print(f"Cubic smoothing spline fitted to {label.upper()} (s = {s_val}).")
except Exception as e:
print(f"{label.upper()} spline failed ({e}), falling back to linear.")
splines[label] = None
spline_fit[label] = False
# Linear fallback interpolators for all three
p50_interp = interp1d(curve_data['wind_bin'], curve_data['p50_power_capped'],
kind='linear', bounds_error=False, fill_value=0.0)
p90_interp = interp1d(curve_data['wind_bin'], curve_data['p90_power_capped'],
kind='linear', bounds_error=False, fill_value=0.0)
p95_interp = interp1d(curve_data['wind_bin'], curve_data['p95_power_capped'],
kind='linear', bounds_error=False, fill_value=0.0)
_fallback = {'p50': p50_interp, 'p90': p90_interp, 'p95': p95_interp}
def power_curve_func(wind_speeds, which='p90'):
"""Return potential power (kW) at given wind speeds.
Physical constraints enforced:
- At or below cut-in (<=2.5 m/s): hard zero
- Normal operating region: spline output capped at 900 kW
- Rated plateau (>=12.5 m/s): forced to 900 kW
- Storm taper (28-34 m/s): linear taper 900 kW to 0 kW
- Above cut-out (>=34 m/s): hard zero
Parameters
----------
wind_speeds : float or array-like
which : 'p50' (lower sensitivity), 'p90' (base case), 'p95' (upper sensitivity)
"""
if which not in ('p50', 'p90', 'p95'):
raise ValueError(f"`which` must be 'p50', 'p90', or 'p95', got {which!r}")
wind_speeds = np.atleast_1d(wind_speeds).astype(float)
if spline_fit[which]:
raw = splines[which](wind_speeds)
else:
raw = _fallback[which](wind_speeds)
result = np.clip(raw, 0, RATED_CAPACITY_KW)
# Cut-in enforcement — hard zero at or below 2.5 m/s
result[wind_speeds <= CUT_IN_MS] = 0.0
# Rated plateau — force 900 kW from rated speed to storm taper start
RATED_SPEED_MS = 12.5
plateau_mask = (wind_speeds > CUT_IN_MS) & \
(wind_speeds >= RATED_SPEED_MS) & \
(wind_speeds < STORM_TAPER_START)
result[plateau_mask] = RATED_CAPACITY_KW
# Storm taper: linear from 900 kW at 28 m/s to 0 kW at 34 m/s
taper_mask = (wind_speeds >= STORM_TAPER_START) & (wind_speeds < CUT_OUT_MS)
result[taper_mask] = RATED_CAPACITY_KW * (
1 - (wind_speeds[taper_mask] - STORM_TAPER_START) /
(CUT_OUT_MS - STORM_TAPER_START)
)
# Full cut-out
result[wind_speeds >= CUT_OUT_MS] = 0.0
return result
# Verify physical constraints at key wind speeds
print(f"\n{'Wind (m/s)':>12} {'P50 (kW)':>10} {'P90 (kW)':>10} {'P95 (kW)':>10}")
for v in [0, 2.0, 2.5, 3, 5, 8, 10, 12, 14, 18, 25, 28, 30, 32, 34, 35]:
p50 = power_curve_func(v, 'p50')[0]
p90 = power_curve_func(v, 'p90')[0]
p95 = power_curve_func(v, 'p95')[0]
print(f" {v:>10.1f} {p50:>10.1f} {p90:>10.1f} {p95:>10.1f}")
Cubic smoothing spline fitted to P50 (s = 3300).
Cubic smoothing spline fitted to P90 (s = 3300).
Cubic smoothing spline fitted to P95 (s = 3300).
Wind (m/s) P50 (kW) P90 (kW) P95 (kW)
0.0 0.0 0.0 0.0
2.0 0.0 0.0 0.0
2.5 0.0 0.0 0.0
3.0 6.0 10.6 12.6
5.0 56.0 76.3 82.7
8.0 254.1 310.6 328.7
10.0 493.4 575.4 600.9
12.0 834.9 892.8 900.0
14.0 900.0 900.0 900.0
18.0 900.0 900.0 900.0
25.0 900.0 900.0 900.0
28.0 900.0 900.0 900.0
30.0 600.0 600.0 600.0
32.0 300.0 300.0 300.0
34.0 0.0 0.0 0.0
35.0 0.0 0.0 0.0
Appendix B.5 — Validation and Sensitivity¶
In [168]:
wind_range = np.linspace(0, 36, 500)
# 5a: Fitted curves over raw scatter — upper envelope check
fig, ax = plt.subplots(figsize=(12, 6))
ax.scatter(uc['Wind_ms'], uc['Power_kw'],
s=0.3, alpha=0.05, c=C5, label='Unconstrained observations')
ax.plot(wind_range, power_curve_func(wind_range, 'p50'),
color=C4, linewidth=1.5, linestyle=':', label='P50 (lower sensitivity)')
ax.plot(wind_range, power_curve_func(wind_range, 'p90'),
color=C1, linewidth=2.5, label='P90 (base case)')
ax.plot(wind_range, power_curve_func(wind_range, 'p95'),
color=C3, linewidth=1.5, linestyle='--', label='P95 (upper sensitivity)')
ax.axhline(y=RATED_CAPACITY_KW, color=C2, linestyle='--',
linewidth=0.8, alpha=0.7, label='Rated 900 kW')
ax.axvline(x=CUT_IN_MS, color=C4, linestyle=':',
linewidth=0.8, alpha=0.7, label=f'Cut-in {CUT_IN_MS} m/s')
ax.axvspan(STORM_TAPER_START, CUT_OUT_MS, color=C2, alpha=0.07)
ax.text(31, 50, 'Storm taper\n28→34 m/s', ha='center', fontsize=9, color=C2)
ax.set_xlabel('Wind speed (m/s)')
ax.set_ylabel('Power output (kW)')
ax.legend(markerscale=6, fontsize=9)
ax.set_xlim(0, 37)
ax.set_ylim(-20, 1000)
# Save appendix version
ax.set_title('Figure B.4: Fitted Power Curves over Unconstrained Scatter')
plt.tight_layout()
plt.savefig('fig_5_5_validation.png', dpi=150, bbox_inches='tight')
# Save main report version
ax.set_title('Figure 5.1: Fitted Power Curves over Unconstrained Scatter')
plt.tight_layout()
plt.savefig('fig_5_1_main_report.png', dpi=150, bbox_inches='tight')
plt.show()
# 5b: Seasonal stability check
df['season'] = df['Timestamp'].dt.month.map({
12: 'Winter', 1: 'Winter', 2: 'Winter',
3: 'Spring', 4: 'Spring', 5: 'Spring',
6: 'Summer', 7: 'Summer', 8: 'Summer',
9: 'Autumn', 10: 'Autumn', 11: 'Autumn'
})
season_colours = {'Winter': C1, 'Spring': C3, 'Summer': C4, 'Autumn': C2}
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(wind_range, power_curve_func(wind_range, 'p90'),
color='black', linewidth=2.5, label='Annual P90 (base case)', zorder=5)
for season, colour in season_colours.items():
s_uc = df[
(df['season'] == season) &
(df['Setpoint_kw'] == RATED_CAPACITY_KW) &
(df['Power_kw'] > 0) &
(df['Wind_ms'] > 0) &
(df['Power_kw'] <= df['Setpoint_kw'])
].copy()
s_uc['wind_bin'] = (s_uc['Wind_ms'] / BIN_WIDTH).round() * BIN_WIDTH
s_pc = s_uc.groupby('wind_bin').agg(
p90=('Power_kw', lambda x: x.quantile(0.90)),
count=('Power_kw', 'count')
).reset_index()
s_pc = s_pc[s_pc['count'] >= MIN_OBS]
s_pc['p90'] = s_pc['p90'].clip(upper=RATED_CAPACITY_KW)
ax.plot(s_pc['wind_bin'], s_pc['p90'],
'o--', color=colour, markersize=3,
linewidth=1.2, alpha=0.8, label=f'{season} P90')
ax.axhline(y=RATED_CAPACITY_KW, color=C2, linestyle='--',
linewidth=0.8, alpha=0.5)
ax.set_xlabel('Wind speed (m/s)')
ax.set_ylabel('Power output (kW)')
ax.set_title('Figure 5.6: Seasonal Stability — P90 by Season vs Annual Curve')
ax.legend(fontsize=9)
ax.set_xlim(0, 30)
ax.set_ylim(-20, 1000)
plt.tight_layout()
plt.savefig('fig_5_6_seasonal_stability.png', dpi=150, bbox_inches='tight')
plt.show()
# 5c: P50/P90/P95 sensitivity — like-for-like using power_curve_func
p50_total_kwh = (power_curve_func(df['Wind_ms'].values, 'p50') * HOURS_PER_STEP).sum()
p90_total_kwh = (power_curve_func(df['Wind_ms'].values, 'p90') * HOURS_PER_STEP).sum()
p95_total_kwh = (power_curve_func(df['Wind_ms'].values, 'p95') * HOURS_PER_STEP).sum()
print(f"\nSensitivity — Potential Energy by Estimator")
print(f" P50 (lower): {p50_total_kwh/1e6:.3f} GWh")
print(f" P90 (base): {p90_total_kwh/1e6:.3f} GWh")
print(f" P95 (upper): {p95_total_kwh/1e6:.3f} GWh")
print(f"\n P90 vs P50: {(p90_total_kwh - p50_total_kwh)/p50_total_kwh*100:+.1f}%")
print(f" P95 vs P90: {(p95_total_kwh - p90_total_kwh)/p90_total_kwh*100:+.1f}%")
--------------------------------------------------------------------------- KeyError Traceback (most recent call last) File /Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/pandas/core/indexes/base.py:3653, in Index.get_loc(self, key) 3652 try: -> 3653 return self._engine.get_loc(casted_key) 3654 except KeyError as err: File /Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/pandas/_libs/index.pyx:147, in pandas._libs.index.IndexEngine.get_loc() File /Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/pandas/_libs/index.pyx:176, in pandas._libs.index.IndexEngine.get_loc() File pandas/_libs/hashtable_class_helper.pxi:7080, in pandas._libs.hashtable.PyObjectHashTable.get_item() File pandas/_libs/hashtable_class_helper.pxi:7088, in pandas._libs.hashtable.PyObjectHashTable.get_item() KeyError: 'Timestamp' The above exception was the direct cause of the following exception: KeyError Traceback (most recent call last) Cell In[168], line 38 35 plt.show() 37 # 5b: Seasonal stability check ---> 38 df['season'] = df['Timestamp'].dt.month.map({ 39 12: 'Winter', 1: 'Winter', 2: 'Winter', 40 3: 'Spring', 4: 'Spring', 5: 'Spring', 41 6: 'Summer', 7: 'Summer', 8: 'Summer', 42 9: 'Autumn', 10: 'Autumn', 11: 'Autumn' 43 }) 45 season_colours = {'Winter': C1, 'Spring': C3, 'Summer': C4, 'Autumn': C2} 47 fig, ax = plt.subplots(figsize=(12, 6)) File /Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/pandas/core/frame.py:3761, in DataFrame.__getitem__(self, key) 3759 if self.columns.nlevels > 1: 3760 return self._getitem_multilevel(key) -> 3761 indexer = self.columns.get_loc(key) 3762 if is_integer(indexer): 3763 indexer = [indexer] File /Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/pandas/core/indexes/base.py:3655, in Index.get_loc(self, key) 3653 return self._engine.get_loc(casted_key) 3654 except KeyError as err: -> 3655 raise KeyError(key) from err 3656 except TypeError: 3657 # If we have a listlike key, _check_indexing_error will raise 3658 # InvalidIndexError. Otherwise we fall through and re-raise 3659 # the TypeError. 3660 self._check_indexing_error(key) KeyError: 'Timestamp'
AppendixPowerCurve.6 — Final Power Curve Outputs¶
In [42]:
# Apply P90 base case to full dataset
df['potential_power_kw'] = power_curve_func(df['Wind_ms'].values, 'p90')
df['potential_energy_kwh'] = df['potential_power_kw'] * HOURS_PER_STEP
df['actual_energy_kwh'] = df['Power_kw'] * HOURS_PER_STEP
df['curtailed_energy_kwh'] = (
(df['potential_power_kw'] - df['Power_kw'])
.clip(lower=0) * HOURS_PER_STEP
)
total_potential = df['potential_energy_kwh'].sum()
total_actual = df['actual_energy_kwh'].sum()
total_curtailed = df['curtailed_energy_kwh'].sum()
total_years = len(df) / (60 * 24 * 365.25)
print("=== Final Potential Power Function — Output Summary ===")
print(f" Data coverage: {total_years:.2f} years ({len(df):,} minutes)")
print(f" Total potential energy: {total_potential/1e6:.3f} GWh")
print(f" Total actual energy: {total_actual/1e6:.3f} GWh")
print(f" Total curtailed energy: {total_curtailed/1e6:.3f} GWh")
print(f" Curtailment fraction: {total_curtailed/total_potential:.1%}")
print(f"\nAnnualised per turbine (P90 base case):")
print(f" Annual potential: {total_potential/total_years/1e3:.0f} MWh/yr")
print(f" Annual actual: {total_actual/total_years/1e3:.0f} MWh/yr")
print(f" Annual curtailed: {total_curtailed/total_years/1e3:.0f} MWh/yr")
# Lookup table at key wind speeds
print(f"\n{'Wind (m/s)':>12} {'P50 (kW)':>10} {'P90 (kW)':>10} {'P95 (kW)':>10}")
for v in [0, 2.5, 4, 6, 8, 10, 12, 14, 18, 24, 28, 30, 32, 34]:
print(f" {v:>10.1f} "
f"{power_curve_func(v, 'p50')[0]:>10.1f} "
f"{power_curve_func(v, 'p90')[0]:>10.1f} "
f"{power_curve_func(v, 'p95')[0]:>10.1f}")
# Export for downstream sections
export_df = df[['Timestamp', 'Power_kw', 'Setpoint_kw', 'Wind_ms', 'state',
'potential_power_kw', 'actual_energy_kwh',
'potential_energy_kwh', 'curtailed_energy_kwh']].copy()
export_df.to_csv('turbine_potential_energy.csv', index=False)
pc_winds = np.arange(0, 40.5, 0.5)
pd.DataFrame({
'wind_speed_ms': pc_winds,
'potential_power_kw_p50': power_curve_func(pc_winds, 'p50'),
'potential_power_kw_p90': power_curve_func(pc_winds, 'p90'),
'potential_power_kw_p95': power_curve_func(pc_winds, 'p95')
}).to_csv('power_curve_lookup.csv', index=False)
print(f"\nExported: turbine_potential_energy.csv ({len(export_df):,} rows)")
print(f"Exported: power_curve_lookup.csv ({len(pc_winds)} rows, P50 + P90 + P95)")
=== Final Potential Power Function — Output Summary ===
Data coverage: 2.02 years (1,064,219 minutes)
Total potential energy: 8.473 GWh
Total actual energy: 6.516 GWh
Total curtailed energy: 1.989 GWh
Curtailment fraction: 23.5%
Annualised per turbine (P90 base case):
Annual potential: 4187 MWh/yr
Annual actual: 3220 MWh/yr
Annual curtailed: 983 MWh/yr
Wind (m/s) P50 (kW) P90 (kW) P95 (kW)
0.0 0.0 0.0 0.0
2.5 0.0 0.0 0.0
4.0 23.9 35.7 40.3
6.0 103.8 134.2 143.0
8.0 254.1 310.6 328.7
10.0 493.4 575.4 600.9
12.0 834.9 892.8 900.0
14.0 900.0 900.0 900.0
18.0 900.0 900.0 900.0
24.0 900.0 900.0 900.0
28.0 900.0 900.0 900.0
30.0 600.0 600.0 600.0
32.0 300.0 300.0 300.0
34.0 0.0 0.0 0.0
Exported: turbine_potential_energy.csv (1,064,219 rows)
Exported: power_curve_lookup.csv (81 rows, P50 + P90 + P95)
Appendix C: Residential Demand Patterns¶
Appendix C.1 — Data Cleaning and Preparation¶
In [169]:
### AppendixDemandPatterns.1 — Data Cleaning and Preparation
import pandas as pd
import numpy as np
# Load raw dataset
df = pd.read_csv('data/Residential_demand.csv')
df['Timestamp'] = pd.to_datetime(df['Timestamp'], format='mixed', dayfirst=True)
df = df.sort_values('Timestamp').reset_index(drop=True)
print(f"Raw dataset: {len(df):,} rows")
print(f"Date range: {df['Timestamp'].min().date()} to {df['Timestamp'].max().date()}")
print(f"\nFirst 5 rows:")
print(df.head())
# Drop single 2018 row
df = df[df['Timestamp'].dt.year == 2017].copy()
df = df.sort_values('Timestamp').reset_index(drop=True)
print(f"\n2017 dataset: {len(df):,} rows (expected {365*48:,})")
# Confirm 30-minute intervals
time_diffs = df['Timestamp'].diff().dropna()
print(f"\nUnique time gaps: {time_diffs.value_counts().to_dict()}")
if len(time_diffs.unique()) == 1:
print("All intervals are exactly 30 minutes. No gap-filling needed.")
# N_households note
print(f"\nN_households range: {df['N_households'].min():,} to {df['N_households'].max():,}")
print("N_households not used in downstream analysis.")
print("Anomaly in Sep-Oct 2017 (~30,000) treated as data quality issue.")
# Outlier check — inspection only, no removal
mean_d = df['Demand_mean_kw'].mean()
std_d = df['Demand_mean_kw'].std()
threshold = mean_d + 3 * std_d
outliers = df[df['Demand_mean_kw'] > threshold]
print(f"\nMean: {mean_d:.4f} kW")
print(f"Std: {std_d:.4f} kW")
print(f"Threshold: {threshold:.4f} kW (mean + 3σ)")
print(f"Outliers flagged: {len(outliers)}")
print(f"Min demand: {df['Demand_mean_kw'].min():.4f} kW")
print(f"Max demand: {df['Demand_mean_kw'].max():.4f} kW")
print(f"Zero or negative values: {(df['Demand_mean_kw'] <= 0).sum()}")
print("Inspection showed flagged observations correspond to expected")
print("winter evening demand peaks and were retained.")
# Aggregate to hourly resolution
df['hour_ts'] = df['Timestamp'].dt.floor('h')
df_hourly = df.groupby('hour_ts').agg(
Demand_mean_kw=('Demand_mean_kw', 'mean')
).reset_index()
df_hourly.rename(columns={'hour_ts': 'Timestamp'}, inplace=True)
energy_30min = (df['Demand_mean_kw'] * 0.5).sum()
energy_hourly = (df_hourly['Demand_mean_kw'] * 1.0).sum()
print(f"\nHourly dataset: {len(df_hourly):,} rows (expected {365*24:,})")
print(f"Energy conservation check:")
print(f" 30-min: {energy_30min:.4f} kWh")
print(f" Hourly: {energy_hourly:.4f} kWh")
print(f" Diff: {abs(energy_30min - energy_hourly):.6f} kWh")
# Assign season and helper columns
# UK meteorological convention:
# Winter=Dec/Jan/Feb, Spring=Mar/Apr/May, Summer=Jun/Jul/Aug, Autumn=Sep/Oct/Nov
def assign_season(month):
if month in [12, 1, 2]: return 'Winter'
elif month in [3, 4, 5]: return 'Spring'
elif month in [6, 7, 8]: return 'Summer'
else: return 'Autumn'
df_hourly['month'] = df_hourly['Timestamp'].dt.month
df_hourly['hour'] = df_hourly['Timestamp'].dt.hour
df_hourly['season'] = df_hourly['month'].apply(assign_season)
print("\nSeason distribution (hourly slots):")
print(df_hourly.groupby('season')['Timestamp'].count())
Raw dataset: 17,568 rows
Date range: 2017-01-01 to 2018-01-01
First 5 rows:
Timestamp Demand_mean_kw N_households
0 2017-01-01 00:00:00 0.220106 5428
1 2017-01-01 00:30:00 0.205945 5429
2 2017-01-01 01:00:00 0.189090 5429
3 2017-01-01 01:30:00 0.173118 5429
4 2017-01-01 02:00:00 0.159051 5429
2017 dataset: 17,520 rows (expected 17,520)
Unique time gaps: {Timedelta('0 days 00:30:00'): 17519}
All intervals are exactly 30 minutes. No gap-filling needed.
N_households range: 5,406 to 30,037
N_households not used in downstream analysis.
Anomaly in Sep-Oct 2017 (~30,000) treated as data quality issue.
Mean: 0.2190 kW
Std: 0.0901 kW
Threshold: 0.4892 kW (mean + 3σ)
Outliers flagged: 21
Min demand: 0.0866 kW
Max demand: 0.5097 kW
Zero or negative values: 0
Inspection showed flagged observations correspond to expected
winter evening demand peaks and were retained.
Hourly dataset: 8,760 rows (expected 8,760)
Energy conservation check:
30-min: 1918.0395 kWh
Hourly: 1918.0395 kWh
Diff: 0.000000 kWh
Season distribution (hourly slots):
season
Autumn 2184
Spring 2208
Summer 2208
Winter 2160
Name: Timestamp, dtype: int64
Appendix C.2 — Seasonal 24-Hour Profiles¶
In [170]:
SEASON_ORDER = ['Winter', 'Spring', 'Summer', 'Autumn']
SEASON_COLORS = {'Winter': C1, 'Spring': C3, 'Summer': C4, 'Autumn': C2}
profiles = (
df_hourly.groupby(['season', 'hour'])['Demand_mean_kw']
.mean()
.reset_index()
)
print("Demand profile summary (kW per household):")
print(profiles.groupby('season')['Demand_mean_kw']
.agg(['min', 'mean', 'max'])
.loc[SEASON_ORDER]
.round(4))
# Peak and trough summary table
peak_summary = []
for season in SEASON_ORDER:
data = profiles[profiles['season'] == season]
peak_row = data.loc[data['Demand_mean_kw'].idxmax()]
trough_row = data.loc[data['Demand_mean_kw'].idxmin()]
peak_summary.append({
'Season' : season,
'Mean (kW)' : round(data['Demand_mean_kw'].mean(), 4),
'Peak time' : f"{int(peak_row['hour']):02d}:00",
'Peak (kW)' : round(peak_row['Demand_mean_kw'], 4),
'Trough time' : f"{int(trough_row['hour']):02d}:00",
'Trough (kW)' : round(trough_row['Demand_mean_kw'], 4)
})
peak_df = pd.DataFrame(peak_summary)
print(peak_df.to_string(index=False))
# Overlaid seasonal profiles — this is the main body figure (Figure 6.1)
tick_hours = list(range(0, 24, 2))
fig, ax = plt.subplots(figsize=(13, 5))
for season in SEASON_ORDER:
data = profiles[profiles['season'] == season].sort_values('hour')
ax.plot(data['hour'], data['Demand_mean_kw'],
color=SEASON_COLORS[season], linewidth=2, label=season)
ax.set_xlabel('Hour of day')
ax.set_ylabel('Average demand per household (kW)')
ax.set_title('Figure 6.1: Seasonal 24-Hour Residential Demand Profiles')
ax.set_xticks(tick_hours)
ax.set_xticklabels([f'{h:02d}:00' for h in tick_hours], rotation=45, fontsize=9)
ax.legend()
plt.tight_layout()
plt.savefig('fig_6_1_seasonal_overlay.png', dpi=150, bbox_inches='tight')
plt.show()
# Four-panel seasonal subplots — appendix figure
fig, axes = plt.subplots(2, 2, figsize=(14, 8), sharey=True)
axes = axes.flatten()
for i, season in enumerate(SEASON_ORDER):
ax = axes[i]
data = profiles[profiles['season'] == season].sort_values('hour')
colour = SEASON_COLORS[season]
ax.plot(data['hour'], data['Demand_mean_kw'],
color=colour, linewidth=2)
ax.fill_between(data['hour'], data['Demand_mean_kw'],
alpha=0.12, color=colour)
peak_row = data.loc[data['Demand_mean_kw'].idxmax()]
trough_row = data.loc[data['Demand_mean_kw'].idxmin()]
ax.axvline(peak_row['hour'], color=colour,
linestyle='--', linewidth=1, alpha=0.7)
ax.annotate(
f"Peak: {int(peak_row['hour']):02d}:00\n{peak_row['Demand_mean_kw']:.3f} kW",
xy=(peak_row['hour'], peak_row['Demand_mean_kw']),
xytext=(peak_row['hour'] - 7, peak_row['Demand_mean_kw'] * 0.96),
fontsize=9, color=colour
)
ax.set_title(season, color=colour)
ax.set_ylabel('Avg demand per household (kW)')
ax.set_xticks(tick_hours)
ax.set_xticklabels([f'{h:02d}:00' for h in tick_hours],
rotation=45, fontsize=9)
fig.suptitle('Figure C.1: Average 24-Hour Demand Profile by Season',
fontsize=12, fontweight='bold')
plt.tight_layout()
plt.savefig('fig_a6_2_seasonal_subplots.png', dpi=150, bbox_inches='tight')
plt.show()
# Heatmap — hour of day vs month
from matplotlib.colors import LinearSegmentedColormap
hso_cmap = LinearSegmentedColormap.from_list('hso', [BG, C4, C1])
heatmap_data = (
df_hourly.groupby(['month', 'hour'])['Demand_mean_kw']
.mean()
.unstack(level='hour')
)
month_labels = ['Jan','Feb','Mar','Apr','May','Jun',
'Jul','Aug','Sep','Oct','Nov','Dec']
fig, ax = plt.subplots(figsize=(14, 5))
im = ax.imshow(heatmap_data.values, aspect='auto',
cmap=hso_cmap, origin='upper')
ax.set_xticks(range(24))
ax.set_xticklabels([f'{h:02d}:00' for h in range(24)],
rotation=45, fontsize=8)
ax.set_yticks(range(12))
ax.set_yticklabels(month_labels)
ax.set_xlabel('Hour of day')
ax.set_ylabel('Month')
ax.set_title('Figure C.2: Average Demand per Household (kW) — Hour vs Month')
plt.colorbar(im, ax=ax, fraction=0.03, pad=0.04,
label='kW per household')
plt.tight_layout()
plt.savefig('fig_a6_3_heatmap.png', dpi=150, bbox_inches='tight')
plt.show()
Demand profile summary (kW per household):
min mean max
season
Winter 0.1128 0.2599 0.4487
Spring 0.0993 0.2115 0.3342
Summer 0.0908 0.1866 0.2845
Autumn 0.0979 0.2187 0.3710
Season Mean (kW) Peak time Peak (kW) Trough time Trough (kW)
Winter 0.2599 18:00 0.4487 04:00 0.1128
Spring 0.2115 18:00 0.3342 03:00 0.0993
Summer 0.1866 17:00 0.2845 03:00 0.0908
Autumn 0.2187 17:00 0.3710 03:00 0.0979
Appendix C.3 — Monthly and Seasonal Demand Variation¶
In [171]:
# Monthly mean per-household demand — shows seasonal shape without scaling
monthly_mean = df_hourly.groupby('month')['Demand_mean_kw'].mean()
month_labels = ['Jan','Feb','Mar','Apr','May','Jun',
'Jul','Aug','Sep','Oct','Nov','Dec']
fig, ax = plt.subplots(figsize=(11, 4))
bars = ax.bar(range(1, 13), monthly_mean.values, color=C1, alpha=0.85)
ax.set_xticks(range(1, 13))
ax.set_xticklabels(month_labels)
ax.set_ylabel('Mean demand per household (kW)')
ax.set_title('Figure C.3: Monthly Mean Residential Demand per Household')
for bar, val in zip(bars, monthly_mean.values):
ax.text(bar.get_x() + bar.get_width() / 2,
bar.get_height() + 0.002,
f'{val:.3f}', ha='center', va='bottom', fontsize=8)
plt.tight_layout()
plt.savefig('fig_a6_4_monthly_demand.png', dpi=150, bbox_inches='tight')
plt.show()
# Seasonal mean demand per household
seasonal_mean = df_hourly.groupby('season')['Demand_mean_kw'].mean()
seasonal_mean = seasonal_mean.loc[SEASON_ORDER]
print("Mean demand per household by season (kW):")
for s, v in seasonal_mean.items():
print(f" {s:<8}: {v:.4f} kW")
print()
print("Note: per-household averages from a representative sample.")
print("Not scaled to total Orkney demand.")
Mean demand per household by season (kW): Winter : 0.2599 kW Spring : 0.2115 kW Summer : 0.1866 kW Autumn : 0.2187 kW Note: per-household averages from a representative sample. Not scaled to total Orkney demand.
Appendix D: Curtailment Analysis¶
Appendix D.1 — Curtailment Calculation Framework¶
In [91]:
from scipy.interpolate import interp1d
# ── Fleet constants ────────────────────────────────────────────────────────────
N_TURBINES = 600
RATED_CAPACITY_KW = 900
EXPORT_CAP_KW = 40_000
# ── Load P90 power curve ───────────────────────────────────────────────────────
pc = pd.read_csv('power_curve_lookup.csv')
pc_interp = interp1d(
pc['wind_speed_ms'].values,
pc['potential_power_kw_p90'].values,
bounds_error=False,
fill_value=0.0
)
print(f'Power curve loaded: {len(pc)} points, '
f'{pc["wind_speed_ms"].min():.1f}–{pc["wind_speed_ms"].max():.1f} m/s')
# ── Load raw 1-minute telemetry ────────────────────────────────────────────────
df_raw = pd.read_csv('data/Turbine_telemetry.csv', parse_dates=['Timestamp'])
df_raw = df_raw.sort_values('Timestamp').reset_index(drop=True)
n_before = len(df_raw)
df_tel = df_raw.dropna(subset=['Power_kw', 'Setpoint_kw', 'Wind_ms']).copy()
print(f'Telemetry: {n_before:,} rows loaded, '
f'{n_before - len(df_tel):,} dropped '
f'({(n_before - len(df_tel)) / n_before * 100:.2f}%)')
print(f'Range: {df_tel["Timestamp"].min()} → {df_tel["Timestamp"].max()}')
# ── Stopped state episode analysis ────────────────────────────────────────────
stopped = df_tel[df_tel['Setpoint_kw'] == 0].copy()
stopped = stopped.sort_values('Timestamp').reset_index(drop=True)
stopped['time_diff'] = stopped['Timestamp'].diff().dt.total_seconds()
stopped['new_episode'] = (stopped['time_diff'] > 120) | (stopped['time_diff'].isna())
stopped['episode_id'] = stopped['new_episode'].cumsum()
episodes = stopped.groupby('episode_id').agg(
duration_hours=('Timestamp', lambda x: (x.max() - x.min()).total_seconds() / 3600),
mean_wind =('Wind_ms', 'mean')
).reset_index()
print(f'\n=== STOPPED EPISODE SUMMARY ===')
print(f' Total stopped episodes: {len(episodes):,}')
print(f'\n {"Duration":<12} {"Episodes":>10} {"Mean wind":>10}')
print(' ' + '-' * 36)
for lo, hi, label in [(0,1,'< 1 hr'), (1,8,'1–8 hr'), (8,24,'8–24 hr'), (24,9999,'> 24 hr')]:
subset = episodes[(episodes['duration_hours'] >= lo) & (episodes['duration_hours'] < hi)]
print(f' {label:<12} {len(subset):>10,} {subset["mean_wind"].mean():>9.1f} m/s')
# ── Apply P90 power curve at 1-minute resolution ───────────────────────────────
df_tel['potential_kw'] = pc_interp(df_tel['Wind_ms'].values)
# ── Resample to hourly means ───────────────────────────────────────────────────
df_curt = (
df_tel.set_index('Timestamp')
[['Power_kw', 'Setpoint_kw', 'Wind_ms', 'potential_kw']]
.resample('1h').mean()
.dropna()
)
print(f'\nHourly series: {len(df_curt):,} rows '
f'({df_curt.index.min()} → {df_curt.index.max()})')
# ── Scale to fleet ─────────────────────────────────────────────────────────────
df_curt['Potential_fleet_kw'] = df_curt['potential_kw'] * N_TURBINES
df_curt['Power_fleet_kw'] = df_curt['Power_kw'] * N_TURBINES
df_curt['Setpoint_fleet_kw'] = df_curt['Setpoint_kw'] * N_TURBINES
# ── Curtailment formula ────────────────────────────────────────────────────────
df_curt['Curtailment_kw'] = (
df_curt['Potential_fleet_kw'] - df_curt['Power_fleet_kw']
).clip(lower=0)
df_curt['Export_kw'] = np.minimum(df_curt['Power_fleet_kw'], EXPORT_CAP_KW)
# ── DR-addressable decomposition ───────────────────────────────────────────────
df_curt['G_kw'] = np.minimum(df_curt['Potential_fleet_kw'],
df_curt['Setpoint_fleet_kw'])
df_curt['Setpoint_loss_kw'] = (df_curt['Potential_fleet_kw']
- df_curt['G_kw']).clip(lower=0)
df_curt['Below_setpoint_kw'] = (df_curt['G_kw']
- df_curt['Power_fleet_kw']).clip(lower=0)
# ── Data span ──────────────────────────────────────────────────────────────────
span_hours_c = (df_curt.index.max() - df_curt.index.min()).total_seconds() / 3600
span_years_c = span_hours_c / 8760
print(f'Data span: {span_years_c:.2f} years ({len(df_curt):,} hourly samples)')
# ── Export ────────────────────────────────────────────────────────────────────
df_curt.to_csv('data/curtailment_hourly.csv')
print(f'Exported: curtailment_hourly.csv ({len(df_curt):,} rows)')
Power curve loaded: 81 points, 0.0–40.0 m/s Telemetry: 1,069,636 rows loaded, 5,417 dropped (0.51%) Range: 2015-05-28 00:00:23 → 2018-01-11 06:14:32 === STOPPED EPISODE SUMMARY === Total stopped episodes: 3,720 Duration Episodes Mean wind ------------------------------------ < 1 hr 3,544 14.4 m/s 1–8 hr 144 14.7 m/s 8–24 hr 24 12.6 m/s > 24 hr 8 10.0 m/s Hourly series: 17,852 rows (2015-05-28 00:00:00 → 2018-01-11 06:00:00) Data span: 2.63 years (17,852 hourly samples) Exported: curtailment_hourly.csv (17,852 rows)
Appendix D.2 — Annual Summary¶
In [98]:
# ── Annualised fleet totals ────────────────────────────────────────────────────
mwh = lambda s: s.sum() / 1000 # hourly kW → MWh
pot_yr = mwh(df_curt['Potential_fleet_kw']) / span_years_c
pwr_yr = mwh(df_curt['Power_fleet_kw']) / span_years_c
exp_yr = mwh(df_curt['Export_kw']) / span_years_c
curt_yr = mwh(df_curt['Curtailment_kw']) / span_years_c
sp_yr = mwh(df_curt['Setpoint_loss_kw']) / span_years_c
bel_yr = mwh(df_curt['Below_setpoint_kw']) / span_years_c
sat_pct = 100 * (df_curt['Power_fleet_kw'] >= EXPORT_CAP_KW - 1).mean()
print('=== ANNUALISED FLEET SUMMARY ===')
print(f' Fleet: {N_TURBINES} × {RATED_CAPACITY_KW} kW = '
f'{N_TURBINES * RATED_CAPACITY_KW / 1000:,.0f} MW installed')
print(f' Span: {span_years_c:.2f} years ({len(df_curt):,} hourly samples)\n')
print(f' Potential: {pot_yr:>12,.0f} MWh/yr ({pot_yr/1000:.1f} GWh/yr)')
print(f' Actual generation: {pwr_yr:>12,.0f} MWh/yr ({100*pwr_yr/pot_yr:.1f}%)')
print(f' Exported (≤40 MW): {exp_yr:>12,.0f} MWh/yr '
f'[info — already embedded in generation]')
print(f' CURTAILED: {curt_yr:>12,.0f} MWh/yr ({100*curt_yr/pot_yr:.1f}%)')
print(f'\n Curtailment value: £{curt_yr * 55:>12,.0f}/yr '
f'@ £55/MWh [generator-side only — not Kaluza revenue]')
print(f'\n DR decomposition:')
print(f' Setpoint-driven (DR-addressable): {sp_yr:>10,.0f} MWh/yr '
f'({100*sp_yr/curt_yr:.1f}%)')
print(f' Below-setpoint (not addressable): {bel_yr:>10,.0f} MWh/yr '
f'({100*bel_yr/curt_yr:.1f}%)')
print(f'\n Cable saturation: {sat_pct:.1f}% of hours at or above 40 MW export')
# ── Monthly stacked bar (Figure 7.1) ──────────────────────────────────────────
# Bar heights are raw observed totals across the full 2.63-year span.
# Months appearing multiple times have proportionally higher bars.
# Partial months at the start and end of the dataset will appear shorter.
# Annualised seasonal comparisons are in AppendixCurtailment.3.
monthly_c = (
df_curt.assign(month=df_curt.index.to_period('M'))
.groupby('month')
.agg(
Generation =('Power_fleet_kw', lambda x: x.sum() / 1e6),
Curtailment =('Curtailment_kw', lambda x: x.sum() / 1e6)
)
)
fig, ax = plt.subplots(figsize=(14, 5))
xp = np.arange(len(monthly_c))
ax.bar(xp, monthly_c['Generation'],
color=C3, label='Actual generation')
ax.bar(xp, monthly_c['Curtailment'],
bottom=monthly_c['Generation'],
color=C2, alpha=0.9, label='Curtailment (stacked height = potential)')
ax.set_xticks(xp)
ax.set_xticklabels([str(p) for p in monthly_c.index],
rotation=45, ha='right', fontsize=8)
ax.set_ylabel('GWh')
ax.set_title(f'Figure 7.1: Monthly Actual Generation and Curtailment '
f'({N_TURBINES}-Turbine Fleet, Observed Data Span)')
ax.legend()
plt.tight_layout()
plt.savefig('fig_7_1_monthly_stacked.png', dpi=150, bbox_inches='tight')
plt.show()
=== ANNUALISED FLEET SUMMARY ===
Fleet: 600 × 900 kW = 540 MW installed
Span: 2.63 years (17,852 hourly samples)
Potential: 1,945,723 MWh/yr (1945.7 GWh/yr)
Actual generation: 1,495,843 MWh/yr (76.9%)
Exported (≤40 MW): 223,989 MWh/yr [info — already embedded in generation]
CURTAILED: 450,431 MWh/yr (23.1%)
Curtailment value: £ 24,773,723/yr @ £55/MWh [generator-side only — not Kaluza revenue]
DR decomposition:
Setpoint-driven (DR-addressable): 268,146 MWh/yr (59.5%)
Below-setpoint (not addressable): 184,271 MWh/yr (40.9%)
Cable saturation: 75.4% of hours at or above 40 MW export
Appendix D.3 — Seasonal Curtailment Patterns¶
In [173]:
def assign_season_c(month):
if month in (12, 1, 2): return 'Winter'
if month in (3, 4, 5): return 'Spring'
if month in (6, 7, 8): return 'Summer'
return 'Autumn'
SEASON_ORDER_C = ['Winter', 'Spring', 'Summer', 'Autumn']
SEASON_COLORS_C = {'Winter': C1, 'Spring': C3, 'Summer': C4, 'Autumn': C2}
SEASON_COLORS_C = {
'Winter': '#243B53', # C5 — dark navy, cold association
'Spring': '#12C2A9', # C3 — teal, fresh
'Summer': '#F15BB5', # Tertiary — warm, expressive
'Autumn': '#7C5CFF' # C4 — secondary, neutral
}
df_curt['Season'] = df_curt.index.month.map(assign_season_c)
df_curt['hour'] = df_curt.index.hour
df_curt['month_num'] = df_curt.index.month
seasonal_c = (
df_curt.groupby('Season')
.agg(
Hours =('Potential_fleet_kw', 'count'),
Potential_MWh =('Potential_fleet_kw', lambda x: x.sum() / 1000),
Power_MWh =('Power_fleet_kw', lambda x: x.sum() / 1000),
Curtailment_MWh =('Curtailment_kw', lambda x: x.sum() / 1000),
Avg_wind_ms =('Wind_ms', 'mean')
)
.reindex(SEASON_ORDER_C)
)
seasonal_c['Potential_MWh'] /= span_years_c
seasonal_c['Power_MWh'] /= span_years_c
seasonal_c['Curtailment_MWh'] /= span_years_c
seasonal_c['Curt_rate_pct'] = (100 * seasonal_c['Curtailment_MWh']
/ seasonal_c['Potential_MWh'])
seasonal_c['Share_annual_pct'] = (100 * seasonal_c['Curtailment_MWh']
/ seasonal_c['Curtailment_MWh'].sum())
print('=== SEASONAL BREAKDOWN (annualised) ===\n')
print(f"{'Season':<8} {'Potential':>14} {'Curtailed':>14} "
f"{'Rate':>8} {'Share':>8} {'Avg wind':>10}")
print('-' * 70)
for s, row in seasonal_c.iterrows():
print(f"{s:<8} {row['Potential_MWh']:>12,.0f} {row['Curtailment_MWh']:>12,.0f} "
f" {row['Curt_rate_pct']:>6.1f}% {row['Share_annual_pct']:>6.1f}% "
f" {row['Avg_wind_ms']:>7.1f} m/s")
# ── Figure 7.2: Seasonal curtailment absolute and rate ────────────────────────
fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))
colors_c = [SEASON_COLORS_C[s] for s in SEASON_ORDER_C]
axes[0].bar(SEASON_ORDER_C, seasonal_c['Curtailment_MWh'] / 1000,
color=colors_c, edgecolor='white', linewidth=0.5)
axes[0].set_ylabel('Curtailment (GWh/yr)')
axes[0].set_title('Figure D.1a: Curtailed Energy by Season (Annualised)')
for i, v in enumerate(seasonal_c['Curtailment_MWh'] / 1000):
axes[0].text(i, v + 0.5, f'{v:,.0f}', ha='center', va='bottom', fontsize=9)
axes[1].bar(SEASON_ORDER_C, seasonal_c['Curt_rate_pct'],
color=colors_c, edgecolor='white', linewidth=0.5)
axes[1].set_ylabel('Curtailment rate (% of potential)')
axes[1].set_title('Figure D.1b: Curtailment Rate by Season')
for i, v in enumerate(seasonal_c['Curt_rate_pct']):
axes[1].text(i, v + 0.3, f'{v:.1f}%', ha='center', va='bottom', fontsize=9)
plt.tight_layout()
plt.savefig('fig_7_2_seasonal.png', dpi=150, bbox_inches='tight')
plt.show()
=== SEASONAL BREAKDOWN (annualised) === Season Potential Curtailed Rate Share Avg wind ---------------------------------------------------------------------- Winter 646,796 76,928 11.9% 17.1% 11.9 m/s Spring 386,608 96,256 24.9% 21.4% 8.9 m/s Summer 372,742 131,930 35.4% 29.3% 8.1 m/s Autumn 539,576 145,317 26.9% 32.3% 10.7 m/s
Appendix D.4 — Daily Curtailment Patterns¶
In [97]:
from matplotlib.colors import LinearSegmentedColormap
# ── Precompute grouped means ───────────────────────────────────────────────────
hourly_mean = df_curt.groupby('hour')['Curtailment_kw'].mean() / 1000
monthly_hourly_mean = df_curt.groupby(['month_num', 'hour'])['Curtailment_kw'].mean() / 1000
# ── Summary statistics ─────────────────────────────────────────────────────────
peak_hour = hourly_mean.idxmax()
peak_mw = hourly_mean.max()
daytime = df_curt[df_curt['hour'].between(8, 17)]['Curtailment_kw'].mean() / 1000
overnight = df_curt[
df_curt['hour'].between(22, 23) | df_curt['hour'].between(0, 6)
]['Curtailment_kw'].mean() / 1000
hot_idx = monthly_hourly_mean.idxmax()
hot_mw = monthly_hourly_mean.max()
month_labels = ['Jan','Feb','Mar','Apr','May','Jun',
'Jul','Aug','Sep','Oct','Nov','Dec']
print('=== DAILY PATTERN SUMMARY ===')
print(f' Peak curtailment hour: {peak_hour:02d}:00 ({peak_mw:.1f} MW mean)')
print(f' Daytime mean (08:00-17:00): {daytime:.1f} MW')
print(f' Overnight mean (22:00-06:00): {overnight:.1f} MW')
print(f' Hottest hour-month cell: '
f'{month_labels[hot_idx[0] - 1]} at {hot_idx[1]:02d}:00 ({hot_mw:.1f} MW)')
# ── Hour × month heatmap (Figure 7.3) ─────────────────────────────────────────
heat = monthly_hourly_mean.unstack('hour')
hso_cmap = LinearSegmentedColormap.from_list('hso', [BG, '#FF5A6E', '#0B1020'])
fig, ax = plt.subplots(figsize=(14, 5))
im = ax.imshow(heat.values, aspect='auto', cmap=hso_cmap,
interpolation='nearest', origin='upper')
ax.set_xticks(range(24))
ax.set_xticklabels([f'{h:02d}:00' for h in range(24)],
rotation=45, fontsize=8)
ax.set_yticks(range(len(heat.index)))
ax.set_yticklabels([month_labels[i - 1] for i in heat.index])
ax.set_xlabel('Hour of day')
ax.set_ylabel('Month')
ax.set_title('Figure 7.2: Mean Curtailment Intensity (MW) '
'by Hour and Month')
plt.colorbar(im, ax=ax, fraction=0.03, pad=0.04, label='Mean curtailment (MW)')
plt.tight_layout()
plt.savefig('fig_7_3_heatmap.png', dpi=150, bbox_inches='tight')
plt.show()
=== DAILY PATTERN SUMMARY === Peak curtailment hour: 13:00 (78.1 MW mean) Daytime mean (08:00-17:00): 71.9 MW Overnight mean (22:00-06:00): 63.6 MW Hottest hour-month cell: Oct at 12:00 (123.5 MW)
In [94]:
# ── Observation counts per month-hour cell ────────────────────────────────────
df_curt['month_num'] = df_curt.index.month
df_curt['hour'] = df_curt.index.hour
month_labels = ['Jan','Feb','Mar','Apr','May','Jun',
'Jul','Aug','Sep','Oct','Nov','Dec']
cell_counts = (df_curt.groupby(['month_num', 'hour'])['Curtailment_kw']
.count()
.unstack('hour'))
obs_per_month = df_curt.groupby('month_num')['Curtailment_kw'].count()
print('=== OBSERVATIONS PER MONTH ===')
print(f' {"Month":<6} {"Hourly rows":>12} {"Avg obs per hour-cell":>22}')
print(' ' + '-' * 44)
for m in range(1, 13):
avg_per_cell = obs_per_month[m] / 24
print(f' {month_labels[m-1]:<6} {obs_per_month[m]:>12,} {avg_per_cell:>22.1f}')
=== OBSERVATIONS PER MONTH === Month Hourly rows Avg obs per hour-cell -------------------------------------------- Jan 1,707 71.1 Feb 1,368 57.0 Mar 1,488 62.0 Apr 949 39.5 May 1,579 65.8 Jun 1,512 63.0 Jul 1,404 58.5 Aug 1,488 62.0 Sep 1,585 66.0 Oct 1,488 62.0 Nov 1,437 59.9 Dec 1,847 77.0
Appendix D.5 — External Plausibility Check¶
In [84]:
ROUSAY_CURT_GWH = 0.7
ROUSAY_SCALED_GWH = ROUSAY_CURT_GWH * N_TURBINES
MODEL_GWH = curt_yr / 1000
ratio = MODEL_GWH / ROUSAY_SCALED_GWH
print('=== SENSE CHECK — ROUSAY TURBINE ===')
print(f' Published Rousay curtailment (FY 2016-17): {ROUSAY_CURT_GWH} GWh')
print(f' Implied fleet total (× {N_TURBINES} turbines): '
f'{ROUSAY_SCALED_GWH:.0f} GWh/yr')
print(f' Model fleet curtailment: '
f'{MODEL_GWH:.0f} GWh/yr')
print(f' Ratio (model / Rousay-implied): {ratio:.2f}×')
=== SENSE CHECK — ROUSAY TURBINE === Published Rousay curtailment (FY 2016-17): 0.7 GWh Implied fleet total (× 600 turbines): 420 GWh/yr Model fleet curtailment: 450 GWh/yr Ratio (model / Rousay-implied): 1.07×
Appendix D.6 — Missing Data Assessment¶
In [174]:
expected_hours = pd.date_range(
df_curt.index.min().floor('h'),
df_curt.index.max().floor('h'),
freq='h'
)
missing_hours = expected_hours.difference(df_curt.index)
print('=== MISSING DATA ASSESSMENT ===')
print(f' Expected hours : {len(expected_hours):,}')
print(f' Present hours : {len(df_curt):,}')
print(f' Missing hours : {len(missing_hours):,} '
f'({100 * len(missing_hours) / len(expected_hours):.1f}%)')
month_labels = ['Jan','Feb','Mar','Apr','May','Jun',
'Jul','Aug','Sep','Oct','Nov','Dec']
exp_by_month = (pd.Series(1, index=expected_hours)
.groupby(expected_hours.month).sum()
.reindex(range(1, 13)))
miss_by_month = (pd.Series(1, index=missing_hours)
.groupby(missing_hours.month).sum()
.reindex(range(1, 13), fill_value=0))
rate_by_month = (miss_by_month / exp_by_month.reindex(range(1, 13)) * 100).round(1)
print(f'\n {"Month":<6} {"Missing hrs":>12} {"Missing rate":>13}')
print(' ' + '-' * 35)
for m in range(1, 13):
print(f' {month_labels[m-1]:<6} {int(miss_by_month[m]):>12,} '
f'{rate_by_month[m]:>12.1f}%')
# ── Missing rate by month bar chart (Figure 7.4) ──────────────────────────────
fig, ax = plt.subplots(figsize=(11, 4))
colors_miss = [C2 if rate_by_month[m] > 25 else C5 for m in range(1, 13)]
ax.bar(range(1, 13), [rate_by_month[m] for m in range(1, 13)],
color=colors_miss, edgecolor='white', linewidth=0.5)
ax.axhline(25, color=C2, linestyle='--', linewidth=1, alpha=0.6,
label='25% threshold')
ax.set_xticks(range(1, 13))
ax.set_xticklabels(month_labels)
ax.set_ylabel('Missing hours (%)')
ax.set_title('Figure D.2: Missing Hourly Data Rate by Month')
from matplotlib.patches import Patch
legend_elements = [
Patch(facecolor=C2, label='Above 25% missing'),
Patch(facecolor=C5, label='Below 25% missing'),
plt.Line2D([0], [0], color=C2, linestyle='--', linewidth=1, alpha=0.6,
label='25% threshold')
]
ax.legend(handles=legend_elements)
plt.tight_layout()
plt.savefig('fig_7_4_missing_data.png', dpi=150, bbox_inches='tight')
plt.show()
=== MISSING DATA ASSESSMENT === Expected hours : 23,023 Present hours : 17,852 Missing hours : 5,171 (22.5%) Month Missing hrs Missing rate ----------------------------------- Jan 28 1.6% Feb 0 0.0% Mar 0 0.0% Apr 491 34.1% May 5 0.3% Jun 648 30.0% Jul 828 37.1% Aug 744 33.3% Sep 575 26.6% Oct 744 33.3% Nov 723 33.5% Dec 385 17.2%
Appendix D.7 — Fleet Size Sensitivity¶
In [175]:
per_turbine_pot_yr = pot_yr / N_TURBINES
per_turbine_pwr_yr = pwr_yr / N_TURBINES
per_turbine_curt_yr = curt_yr / N_TURBINES
fleet_sizes = [500, 600, 700]
print('=== FLEET SIZE SENSITIVITY ===')
print(f' Per-turbine annual potential: {per_turbine_pot_yr:,.0f} MWh/yr')
print(f' Per-turbine annual curtailment: {per_turbine_curt_yr:,.0f} MWh/yr')
print(f' Per-turbine curtailment rate: {100*per_turbine_curt_yr/per_turbine_pot_yr:.1f}%')
print()
print(f' {"Fleet":>8} {"Installed MW":>14} {"Potential GWh/yr":>18} '
f'{"Curtailed GWh/yr":>18} {"Curt rate":>10}')
print(' ' + '-' * 75)
for n in fleet_sizes:
pot_s = per_turbine_pot_yr * n / 1000
curt_s = per_turbine_curt_yr * n / 1000
rate_s = 100 * curt_s / pot_s
marker = ' ← base case' if n == N_TURBINES else ''
print(f' {n:>8} {n * RATED_CAPACITY_KW / 1000:>14,.0f} '
f'{pot_s:>18,.0f} {curt_s:>18,.0f} {rate_s:>9.1f}%{marker}')
# ── Fleet sensitivity bar chart (Figure 7.5) ──────────────────────────────────
curt_vals = [per_turbine_curt_yr * n / 1000 for n in fleet_sizes]
bar_colors = [C5, C3, C4]
fig, ax = plt.subplots(figsize=(8, 4))
bars = ax.bar([str(n) for n in fleet_sizes], curt_vals,
color=bar_colors, edgecolor='white', linewidth=0.5)
ax.set_xlabel('Number of turbines')
ax.set_ylabel('Annual curtailment (GWh/yr)')
ax.set_title('Figure D.3: Fleet Size Sensitivity of Annual Curtailment Estimates')
for bar, val in zip(bars, curt_vals):
ax.text(bar.get_x() + bar.get_width() / 2,
val + 2, f'{val:,.0f}',
ha='center', va='bottom', fontsize=10)
plt.tight_layout()
plt.savefig('fig_7_5_fleet_sensitivity.png', dpi=150, bbox_inches='tight')
plt.show()
=== FLEET SIZE SENSITIVITY ===
Per-turbine annual potential: 3,243 MWh/yr
Per-turbine annual curtailment: 751 MWh/yr
Per-turbine curtailment rate: 23.1%
Fleet Installed MW Potential GWh/yr Curtailed GWh/yr Curt rate
---------------------------------------------------------------------------
500 450 1,621 375 23.1%
600 540 1,946 450 23.1% ← base case
700 630 2,270 526 23.1%
Appendix F: Financial Analysis¶
Appendix F.1 — Cost Structure¶
In [130]:
# ============================================================
# AppendixFinancialModel.1 — Cost Structure
# Corresponds to Section 9.1 of the main report.
# Assumes HSO style constants (C1-C5, BG) already defined.
# ============================================================
import numpy as np
import pandas as pd
# ── Market parameters ────────────────────────────────────────────────────────
N_HH = 10_385 # Total Orkney households
# Source: Project FAQ
SCA_SHARE = 0.32 # Sc.A — existing storage heater households
SCB_SHARE = 0.28 # Sc.B — oil/coal households
# Source: OREF EPC survey 2022 — used as 2017/18 proxy
# https://www.oref.co.uk/orkneys-energy/energy-in-orkney/
FUEL_POV = 0.58 # Orkney-specific fuel poverty rate
# Source: Orkney Sustainable Energy Strategy 2016
# https://www.orkney.gov.uk/media/ga4fqvno/sustainable_orkney_energy_strategy_accessible.pdf
N_SCA = int(np.floor(N_HH * SCA_SHARE + 0.5)) # 3,323
N_SCB = int(np.floor(N_HH * SCB_SHARE + 0.5)) # 2,908
# ── Device and installation costs ────────────────────────────────────────────
DEVICE_A = 100 # £/HH — retrofit smart hub
# Source: Project FAQ; Energy Brief 2026
DEVICE_B = 1_200 # £/HH — 2× Dimplex Quantum (2 × £600)
# Source: Energy Brief 2026; HeaterShop 2017
INSTALL_A = 125 # £/HH base — ~17% downward adjustment from Refurbb guide
# https://www.refurbb.co.uk
# Reflects 2017/18 labour market conditions
INSTALL_A_HI = 150 # £/HH upper bound sensitivity
INSTALL_B = 250 # £/HH — full replacement, 2 heaters + wiring
# Internal modelling; Checkatrade current adjusted
# https://www.checkatrade.com
# No pre-2017 primary source — flagged in Section 12
CAPEX_A = DEVICE_A + INSTALL_A # £225/HH base
CAPEX_B = DEVICE_B + INSTALL_B # £1,450/HH
DEVICE_LIFE = 12 # years — Dimplex Quantum 10-15yr manufacturer spec
# https://www.dimplex.co.uk
# ── Grant structure ───────────────────────────────────────────────────────────
# Warmer Homes Scotland — covers full device + install for eligible HH
# Average grant £4,572/HH
# Source: Scottish Govt HEEP Annual Review 2016/17
# https://www.gov.scot/publications/home-energy-efficiency-programmes-scotland-warmer-homes-scotland-2016-17/pages/4/
# Eligibility based on Orkney fuel poverty rate of 58%
# Additional subsidy scenarios (0%, 25%, 50%, 75%) evaluated in Section 10
N_SCA_ELIG = int(N_SCA * FUEL_POV)
N_SCA_NONELIG = N_SCA - N_SCA_ELIG
N_SCB_ELIG = int(N_SCB * FUEL_POV)
N_SCB_NONELIG = N_SCB - N_SCB_ELIG
# ── Operating costs ───────────────────────────────────────────────────────────
MAINT = 31 # £/HH/yr — internal modelling
# No published benchmark identified
# Conservative relative to smart home service plans (~£99/yr)
# Flagged in Section 12
KALUZA_INC = 48 # £/HH/yr — Kaluza portion (60% of £80 base)
GOVT_INC = 32 # £/HH/yr — Government portion (40% of £80 base)
HH_BENEFIT = KALUZA_INC + GOVT_INC # £80/HH/yr total base case
HH_BENEFIT_HI = 100 # £/HH/yr upper bound sensitivity
# Note: no 2017 comparator available — internal modelling
# Flagged in Section 12 as key uncertainty
OPEX_PER_HH = MAINT + KALUZA_INC # £79/HH/yr Kaluza bears
ACQ = 40 # £/HH — customer acquisition, internal modelling
ONBOARD = HH_BENEFIT * 3 / 12 # £20/HH — 3 months free incentive
ONE_OFF = ACQ + ONBOARD # £60/HH total one-off
# ── Fixed annual overhead ─────────────────────────────────────────────────────
# Salary benchmarks: ONS ASHE 2017
# https://www.ons.gov.uk/employmentandlabourmarket/peopleinwork/earningsandworkinghours/bulletins/annualsurveyofhoursandearnings/2017provisionaland2016revisedresults
#
# Employer NI: 13.8% on earnings above £157/week secondary threshold
# Source: HMRC 2017/18
# https://www.gov.uk/employer-national-insurance-contributions
#
# Employer pension: 3% auto-enrolment minimum 2017
# https://www.gov.uk/workplace-pensions/what-you-your-employer-and-the-government-pay
OPS_SALARY = 35_000 # £/yr — below ONS ASHE median
# Reflects small organisation and remote Orkney location
OPS_NI = OPS_SALARY * 0.138
OPS_PENSION = OPS_SALARY * 0.03
OPS_LOADED = OPS_SALARY + OPS_NI + OPS_PENSION # £40,880/yr
SUP_SALARY = 25_000 # £/yr — above customer service median
# Reflects technical role requirements
SUP_NI = SUP_SALARY * 0.138
SUP_PENSION = SUP_SALARY * 0.03
SUP_LOADED = SUP_SALARY + SUP_NI + SUP_PENSION # £29,600/yr
PLATFORM = 7_500 # £/yr — platform hosting and monitoring, internal modelling
TRAVEL = 7_520 # £/yr — travel, admin, local liaison
# Internal modelling — reflects Orkney island logistics
FIXED_OPEX = round(OPS_LOADED + SUP_LOADED + PLATFORM + TRAVEL)
# ── Output tables ─────────────────────────────────────────────────────────────
print("AppendixFinancialModel.1 — Cost Structure")
print("=" * 72)
print(f"\n── Device and Installation Costs ──")
print(f"\n {'Item':<42} {'Value':>12} Source")
print(" " + "-" * 72)
rows = [
('Sc.A device (smart hub)', f'£{DEVICE_A}/HH', 'Project FAQ; Energy Brief 2026'),
('Sc.B device (2× Dimplex Quantum)', f'£{DEVICE_B:,}/HH', 'Energy Brief 2026; HeaterShop 2017'),
('Sc.A installation (base)', f'£{INSTALL_A}/HH', 'Refurbb guide, 2017-adjusted'),
('Sc.A installation (upper bound)', f'£{INSTALL_A_HI}/HH', 'Sensitivity'),
('Sc.B installation', f'£{INSTALL_B}/HH', 'Internal modelling; Checkatrade adjusted'),
('Sc.A total capex (base)', f'£{CAPEX_A}/HH', 'Derived'),
('Sc.B total capex', f'£{CAPEX_B:,}/HH', 'Derived'),
('Device lifetime', f'{DEVICE_LIFE} years', 'Dimplex Quantum manufacturer spec'),
]
for item, val, src in rows:
print(f" {item:<42} {val:>12} {src}")
print(f"\n── Grant Structure ──")
print(f" Fuel poverty rate (Orkney): {FUEL_POV*100:.0f}%")
print(f" Source: Orkney Sustainable Energy Strategy 2016")
print(f" Warmer Homes Scotland avg grant: £4,572/HH")
print(f" Source: Scottish Govt HEEP Annual Review 2016/17")
print(f" Sc.A eligible (fully grant-funded): ~{N_SCA_ELIG:,} HH")
print(f" Sc.A non-eligible (Kaluza-funded): ~{N_SCA_NONELIG:,} HH × £{CAPEX_A} = £{N_SCA_NONELIG*CAPEX_A:,}")
print(f" Sc.B eligible (fully grant-funded): ~{N_SCB_ELIG:,} HH")
print(f" Sc.B non-eligible (Kaluza-funded): ~{N_SCB_NONELIG:,} HH × £{CAPEX_B:,} = £{N_SCB_NONELIG*CAPEX_B:,}")
print(f" Additional subsidy scenarios (0-75%): Evaluated in Section 10")
print(f"\n── Operating Costs ──")
print(f"\n {'Item':<42} {'Value':>12} Source")
print(" " + "-" * 72)
opex_rows = [
('Maintenance', f'£{MAINT}/HH/yr', 'Internal modelling — Section 12'),
('Participation incentive (base)', f'£{HH_BENEFIT}/HH/yr', 'Internal modelling — Section 12'),
(' Kaluza portion (60%)', f'£{KALUZA_INC}/HH/yr', 'Internal modelling'),
(' Government portion (40%)', f'£{GOVT_INC}/HH/yr', 'Internal modelling'),
('Incentive upper bound', f'£{HH_BENEFIT_HI}/HH/yr','Sensitivity'),
('Kaluza opex per HH', f'£{OPEX_PER_HH}/HH/yr', 'Derived'),
('Customer acquisition', f'£{ACQ}/HH', 'Internal modelling'),
('Onboarding (3 months free)', f'£{ONBOARD:.0f}/HH', 'Derived: 3 × (£80/12)'),
('Total one-off per HH', f'£{ONE_OFF:.0f}/HH', 'Derived'),
]
for item, val, src in opex_rows:
print(f" {item:<42} {val:>12} {src}")
print(f"\n── Fixed Annual Overhead: £{FIXED_OPEX:,}/yr ──")
print(f"\n {'Item':<30} {'Salary':>10} {'NI':>8} {'Pension':>8} {'Total':>10} Source")
print(" " + "-" * 76)
staff = [
('Operations manager', OPS_SALARY, OPS_NI, OPS_PENSION, OPS_LOADED, 'ONS ASHE 2017'),
('Customer support', SUP_SALARY, SUP_NI, SUP_PENSION, SUP_LOADED, 'ONS ASHE 2017'),
]
for name, sal, ni, pen, total, src in staff:
print(f" {name:<30} £{sal:>8,} £{ni:>6,.0f} £{pen:>6,.0f} £{total:>8,.0f} {src}")
print(f" {'Platform hosting':<30} {'':>10} {'':>8} {'':>8} £{PLATFORM:>8,} Internal modelling")
print(f" {'Travel, admin, liaison':<30} {'':>10} {'':>8} {'':>8} £{TRAVEL:>8,} Internal modelling")
print(f" {'TOTAL':<30} {'':>10} {'':>8} {'':>8} £{FIXED_OPEX:>8,}")
print(f"\n── Island Logistics Note ──")
print(f" Orkney's remote archipelago geography introduces ferry delivery,")
print(f" limited local contractor availability, and elevated travel costs.")
print(f" Installation assumptions are conservative estimates.")
print(f" Full Sc.A sensitivity (£100-£175/HH) evaluated in Section 10.")
AppendixFinancialModel.1 — Cost Structure
========================================================================
── Device and Installation Costs ──
Item Value Source
------------------------------------------------------------------------
Sc.A device (smart hub) £100/HH Project FAQ; Energy Brief 2026
Sc.B device (2× Dimplex Quantum) £1,200/HH Energy Brief 2026; HeaterShop 2017
Sc.A installation (base) £125/HH Refurbb guide, 2017-adjusted
Sc.A installation (upper bound) £150/HH Sensitivity
Sc.B installation £250/HH Internal modelling; Checkatrade adjusted
Sc.A total capex (base) £225/HH Derived
Sc.B total capex £1,450/HH Derived
Device lifetime 12 years Dimplex Quantum manufacturer spec
── Grant Structure ──
Fuel poverty rate (Orkney): 58%
Source: Orkney Sustainable Energy Strategy 2016
Warmer Homes Scotland avg grant: £4,572/HH
Source: Scottish Govt HEEP Annual Review 2016/17
Sc.A eligible (fully grant-funded): ~1,927 HH
Sc.A non-eligible (Kaluza-funded): ~1,396 HH × £225 = £314,100
Sc.B eligible (fully grant-funded): ~1,686 HH
Sc.B non-eligible (Kaluza-funded): ~1,222 HH × £1,450 = £1,771,900
Additional subsidy scenarios (0-75%): Evaluated in Section 10
── Operating Costs ──
Item Value Source
------------------------------------------------------------------------
Maintenance £31/HH/yr Internal modelling — Section 12
Participation incentive (base) £80/HH/yr Internal modelling — Section 12
Kaluza portion (60%) £48/HH/yr Internal modelling
Government portion (40%) £32/HH/yr Internal modelling
Incentive upper bound £100/HH/yr Sensitivity
Kaluza opex per HH £79/HH/yr Derived
Customer acquisition £40/HH Internal modelling
Onboarding (3 months free) £20/HH Derived: 3 × (£80/12)
Total one-off per HH £60/HH Derived
── Fixed Annual Overhead: £85,100/yr ──
Item Salary NI Pension Total Source
----------------------------------------------------------------------------
Operations manager £ 35,000 £ 4,830 £ 1,050 £ 40,880 ONS ASHE 2017
Customer support £ 25,000 £ 3,450 £ 750 £ 29,200 ONS ASHE 2017
Platform hosting £ 7,500 Internal modelling
Travel, admin, liaison £ 7,520 Internal modelling
TOTAL £ 85,100
── Island Logistics Note ──
Orkney's remote archipelago geography introduces ferry delivery,
limited local contractor availability, and elevated travel costs.
Installation assumptions are conservative estimates.
Full Sc.A sensitivity (£100-£175/HH) evaluated in Section 10.
Appendix F.2 — Revenue Streams¶
In [176]:
# ============================================================
# AppendixFinancialModel.2 — Revenue Streams
# Corresponds to Section 9.2 of the main report.
# Assumes HSO style constants (C1-C5, BG) and cost constants
# from AppendixFinancialModel.1 already defined.
# ============================================================
# ── Revenue constants ────────────────────────────────────────────────────────
SSEN_BASE = 250 # £/MWh base case — 17% discount from £300 mature benchmark
# Reflects early-stage CMZ market in 2017/18
# Source: TED notice 292933-2017 (27 July 2017)
# https://ted.europa.eu/en/notice/-/detail/292933-2017
# SSEN CMZ procurement launch 2016
# https://www.ssen.co.uk/news-views/2016/2016-ssen-opens-constraint-managed-zone/
SSEN_UPSIDE = 300 # £/MWh upside sensitivity — mature market benchmark
# Source: The Energyst, March 2019
# https://theenergyst.com/ssen-procure-flex-across-entire-network-households-evs/
SSEN_SCENS = [150, 200, 250, 300, 400]
GEN_FEE_BASE = 25 # £/MWh — generator service fee base case, internal modelling
GEN_FEE_SCENS = [0, 15, 25, 40]
WHOLESALE = 50 # £/MWh base — BEIS QEP H1 2018
# https://www.gov.uk/government/collections/quarterly-energy-prices
WHOLESALE_LO = 45 # £/MWh downside sensitivity
WHOLESALE_HI = 55 # £/MWh upside sensitivity
ROC_VALUE = 41 # £/MWh — 0.9 ROCs/MWh × £45.58 Ofgem buy-out price 2017/18
# Source: Ofgem RO Guidance — Suppliers
# https://www.ofgem.gov.uk/sites/default/files/2025-03/Renewables-Obligation-(RO)-Guidance-Suppliers.pdf
E7_OFF = 7.5 / 100 # £/kWh — derived from Ofgem prepayment cap Apr-Sep 2018
# See derivation in AppendixFinancialModel.3
# ── Revenue summary table ─────────────────────────────────────────────────────
print("AppendixFinancialModel.2 — Revenue Streams")
print("=" * 72)
print(f"\n── Primary Revenue: SSEN-Style Flexibility Payment ──")
print(f" Base case: £{SSEN_BASE}/MWh [TED notice 292933-2017; SSEN CMZ launch 2016]")
print(f" Upside: £{SSEN_UPSIDE}/MWh [The Energyst, March 2019 — mature market benchmark]")
print(f" Sensitivity: £{min(SSEN_SCENS)}-£{max(SSEN_SCENS)}/MWh")
print(f" Note: Benchmark proxy — not a guaranteed contract rate")
print(f" Rationale: 17% discount from £300 reflects early-stage CMZ market")
print(f" in 2017/18 when provider competition was limited")
print(f"\n── Secondary Revenue: Generator Service Fee ──")
print(f" Base case: £{GEN_FEE_BASE}/MWh [Internal modelling]")
print(f" Range: £{min(GEN_FEE_SCENS)}-£{max(GEN_FEE_SCENS)}/MWh")
print(f"\n── Household Participation Incentive ──")
print(f" Base case: £{HH_BENEFIT}/HH/yr (Kaluza £{KALUZA_INC} + Govt £{GOVT_INC})")
print(f" Upper bound: £{HH_BENEFIT_HI}/HH/yr")
print(f" Note: Funded from institutional revenues — no household subscription")
# ── Generator value check table ──────────────────────────────────────────────
print(f"\n── Generator Value Check (£/MWh of unlocked curtailed generation) ──")
print(f"\n {'Scenario':<28} {'Wholesale':>10} {'ROC':>8} {'Total':>8} "
f"{'Fee':>8} {'Net':>8} {'Rational?':>12}")
print(" " + "-" * 88)
scenarios = [
('Base case (wholesale only)', WHOLESALE, 0, GEN_FEE_BASE),
('Wholesale lower', WHOLESALE_LO, 0, GEN_FEE_BASE),
('Wholesale upper', WHOLESALE_HI, 0, GEN_FEE_BASE),
('ROC upside (accredited)', WHOLESALE, ROC_VALUE, GEN_FEE_BASE),
('Max fee', WHOLESALE, 0, max(GEN_FEE_SCENS)),
]
for label, ws, roc, fee in scenarios:
total = ws + roc
net = total - fee
flag = '✓ rational' if ws - fee > 0 else '✗ irrational'
print(f" {label:<26} £{ws:>6}/MWh £{roc:>4}/MWh £{total:>4}/MWh "
f"£{fee:>4}/MWh £{net:>4}/MWh {flag}")
print(f"\n Note: ROC-inclusive case applies to turbines accredited under")
print(f" the Renewables Obligation before 2015. Base case uses wholesale")
print(f" only (net £{WHOLESALE - GEN_FEE_BASE}/MWh). Curtailed energy currently earns £0/MWh.")
print(f" Participation is commercially rational in all scenarios tested.")
# ── Generator value chart ─────────────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(10, 5))
categories = ['Curtailed\n(baseline)', 'Wholesale\nrevenue',
'Service fee\n(paid to Kaluza)', 'Net to\nwind farmer']
values = [0, WHOLESALE, GEN_FEE_BASE, WHOLESALE - GEN_FEE_BASE]
bottoms = [0, 0, WHOLESALE - GEN_FEE_BASE, 0]
colours = [C4, C3, C2, C1]
labels = ['£0', f'£{WHOLESALE}', f'-£{GEN_FEE_BASE}',
f'£{WHOLESALE - GEN_FEE_BASE}']
for i in range(4):
if i == 2:
ax.bar(i, values[i], bottom=bottoms[i],
color=colours[i], alpha=0.85, width=0.5)
ax.text(i, bottoms[i] - 4, labels[i],
ha='center', fontsize=9, fontweight='bold', color=C2)
ax.text(i, bottoms[i] + values[i] / 2,
'paid to\nKaluza', ha='center', fontsize=7,
color='white', va='center')
elif i == 3:
ax.bar(i, values[i], bottom=0,
color=colours[i], alpha=0.85, width=0.5)
ax.text(i, values[i] + 1.5, labels[i],
ha='center', fontsize=9, fontweight='bold', color='#111111')
else:
ax.bar(i, values[i], bottom=bottoms[i],
color=colours[i], alpha=0.85, width=0.5)
ax.text(i, bottoms[i] + values[i] + 1.5, labels[i],
ha='center', fontsize=9, fontweight='bold', color='#111111')
# ROC upside annotation
ax.annotate(f'ROC-accredited turbines\nretain £{WHOLESALE + ROC_VALUE - GEN_FEE_BASE}/MWh (upside)',
xy=(3, WHOLESALE - GEN_FEE_BASE),
xytext=(3.3, WHOLESALE - GEN_FEE_BASE + 20),
fontsize=8, color='#444444',
arrowprops=dict(arrowstyle='->', color='#444444', lw=0.8))
ax.axhline(0, color='#444444', linewidth=0.6)
ax.set_xticks(range(4))
ax.set_xticklabels(categories, fontsize=10)
ax.set_ylabel('£ per MWh of unlocked generation')
ax.set_title(
'Figure 9.1: Wind Farmer Value per MWh of Unlocked Curtailed Generation',
fontsize=11,
fontweight='bold',
pad=20
)
ax.text(
0.5, 1.02,
'Base case: wholesale £50/MWh, service fee £25/MWh',
transform=ax.transAxes,
ha='center',
va='bottom',
fontsize=9,
color='dimgray'
)
ax.set_ylim(-15, 75)
plt.tight_layout()
plt.savefig('fig_9_1_generator_value.png', dpi=150, bbox_inches='tight')
plt.show()
AppendixFinancialModel.2 — Revenue Streams
========================================================================
── Primary Revenue: SSEN-Style Flexibility Payment ──
Base case: £250/MWh [TED notice 292933-2017; SSEN CMZ launch 2016]
Upside: £300/MWh [The Energyst, March 2019 — mature market benchmark]
Sensitivity: £150-£400/MWh
Note: Benchmark proxy — not a guaranteed contract rate
Rationale: 17% discount from £300 reflects early-stage CMZ market
in 2017/18 when provider competition was limited
── Secondary Revenue: Generator Service Fee ──
Base case: £25/MWh [Internal modelling]
Range: £0-£40/MWh
── Household Participation Incentive ──
Base case: £80/HH/yr (Kaluza £48 + Govt £32)
Upper bound: £100/HH/yr
Note: Funded from institutional revenues — no household subscription
── Generator Value Check (£/MWh of unlocked curtailed generation) ──
Scenario Wholesale ROC Total Fee Net Rational?
----------------------------------------------------------------------------------------
Base case (wholesale only) £ 50/MWh £ 0/MWh £ 50/MWh £ 25/MWh £ 25/MWh ✓ rational
Wholesale lower £ 45/MWh £ 0/MWh £ 45/MWh £ 25/MWh £ 20/MWh ✓ rational
Wholesale upper £ 55/MWh £ 0/MWh £ 55/MWh £ 25/MWh £ 30/MWh ✓ rational
ROC upside (accredited) £ 50/MWh £ 41/MWh £ 91/MWh £ 25/MWh £ 66/MWh ✓ rational
Max fee £ 50/MWh £ 0/MWh £ 50/MWh £ 40/MWh £ 10/MWh ✓ rational
Note: ROC-inclusive case applies to turbines accredited under
the Renewables Obligation before 2015. Base case uses wholesale
only (net £25/MWh). Curtailed energy currently earns £0/MWh.
Participation is commercially rational in all scenarios tested.
Appendix F.3 — Household Economics¶
In [ ]:
# ============================================================
#
# Corresponds to Section 9.3 of the main report.
# Assumes HSO style constants (C1-C5, BG) and constants
# from AppendixFinancialModel.1 and .2 already defined.
# ============================================================
# ── Economy 7 tariff derivation ───────────────────────────────────────────────
# Source: Ofgem prepayment price cap Apr-Sep 2018, North Scotland Economy 7
# https://www.ofgem.gov.uk/system/files/docs/2018/02/prepayment_price_cap_-_01_april_2018_to_30_september_2018_0.pdf
#
# Total annual bill (4,600 kWh): £687.14
# Standing charge: £111.38/yr
# Energy cost: £575.76
#
# Day/night split from Ofgem TDCV August 2017 (Profile Class 2):
# Day: 4,600 × 0.58 = 2,668 kWh
# Night: 4,600 × 0.42 = 1,932 kWh
#
# Assumed day rate ~16p/kWh (standard E7 structure):
# £575.76 = (2,668 × 0.16) + (1,932 × night rate)
# £575.76 = £426.88 + (1,932 × night rate)
# Night rate = £148.88 / 1,932 = 7.7p/kWh → rounded to 7.5p/kWh
#
# FLAGGED: derived figure, not directly stated. Treat as approximate.
# Section 12 limitation.
E7_PEAK = 16.0 / 100 # £/kWh — assumed day rate (standard E7 structure)
E7_OFF = 7.5 / 100 # £/kWh — derived night rate (see derivation above)
# ── Oil price and efficiency ──────────────────────────────────────────────────
OIL_P_L = 39.17 / 100 # £/litre — BEIS QEP Table 4.1.2, 2017 annual average
# https://www.gov.uk/government/collections/quarterly-energy-prices
CV_OIL = 10.18 # kWh/litre — kerosene calorific value
# Source: BEIS greenhouse gas conversion factors 2017
# https://www.gov.uk/government/publications/greenhouse-gas-reporting-conversion-factors-2017
BOIL_EFF = 0.85 # oil boiler efficiency — HHIC industry standard
# https://www.hhic.org.uk
OIL_EFF = OIL_P_L / (CV_OIL * BOIL_EFF) # £/kWh useful heat = 4.527p/kWh
# ── Sc.A consumption ─────────────────────────────────────────────────────────
# Ofgem TDCV August 2017, Profile Class 2 (Economy 7 medium)
# https://www.ofgem.gov.uk/system/files/docs/2017/08/tdcvs_2017_decision.pdf
# Base: 4,200 kWh/yr × 1.2 Orkney climate uplift = 5,040 kWh/yr
# Uplift: internal modelling — no published source, flagged in Section 12
# Day/night split: 58% day / 42% night per Ofgem TDCV 2017
SCA_KWH = 4200 * 1.20
SCA_NIGHT = SCA_KWH * 0.42
SCA_DAY = SCA_KWH * 0.58
SCA_COST = SCA_NIGHT * E7_OFF + SCA_DAY * E7_PEAK
# ── Sc.B consumption ─────────────────────────────────────────────────────────
# Annual heating demand: 12,000 kWh/yr
# Proxy: Ofgem TDCV August 2017 medium gas consumption value (12,000 kWh/yr)
# applied as indicative total heating demand for a medium off-gas Scottish home
# No specific EST 2017 source confirmed — flagged in Section 12
SCB_KWH = 12_000
SCB_OIL_COST = SCB_KWH * OIL_EFF
SCB_ELEC = SCB_KWH * E7_OFF
SCB_NET = SCB_OIL_COST - SCB_ELEC + HH_BENEFIT
# ── Summary table ────────────────────────────────────────────────────────────
print("AppendixFinancialModel.3 — Household Economics")
print("=" * 72)
print(f"\n── Sc.A — Existing Storage Heater Household ──")
print(f" Annual consumption: {SCA_KWH:.0f} kWh/yr")
print(f" Ofgem TDCV Aug 2017 (4,200 kWh) × 1.2 Orkney uplift")
print(f" Night consumption: {SCA_NIGHT:.0f} kWh/yr (42%)")
print(f" Day consumption: {SCA_DAY:.0f} kWh/yr (58%)")
print(f" Night cost: £{SCA_NIGHT * E7_OFF:.2f} ({E7_OFF*100}p/kWh E7 off-peak, derived)")
print(f" Day cost: £{SCA_DAY * E7_PEAK:.2f} ({E7_PEAK*100}p/kWh E7 peak, derived)")
print(f" Current annual cost: £{SCA_COST:.2f}/yr (£{SCA_COST/12:.2f}/month)")
print(f" Participation incentive: +£{HH_BENEFIT}/yr")
print(f" Upfront device cost: £0 (Kaluza-funded or grant-covered)")
print(f" Net household benefit: +£{HH_BENEFIT}/yr plus potential scheduling saving")
print(f"\n── Sc.B — Oil/Coal Household Converting to Storage Heaters ──")
print(f" Annual heating demand: {SCB_KWH:,} kWh/yr")
print(f" Ofgem TDCV Aug 2017 gas medium proxy — flagged in Section 12")
print(f" Oil price (2017 avg): {OIL_P_L*100:.2f}p/litre (BEIS QEP Table 4.1.2, 2017)")
print(f" Calorific value: {CV_OIL} kWh/litre (BEIS conversion factors 2017)")
print(f" Boiler efficiency: {BOIL_EFF*100:.0f}% (HHIC industry standard)")
print(f" Oil effective cost: {OIL_EFF*100:.3f}p/kWh useful heat (derived)")
print(f" Current oil cost: £{SCB_OIL_COST:.2f}/yr (£{SCB_OIL_COST/12:.2f}/month)")
print(f" New E7 electricity: £{SCB_ELEC:.2f}/yr ({SCB_KWH:,} kWh × {E7_OFF*100}p)")
print(f" Participation incentive: +£{HH_BENEFIT}/yr")
print(f" Net vs oil baseline: £{SCB_NET:.2f}/yr ({'saving' if SCB_NET > 0 else 'additional cost'})")
print(f"\n── Household Economics Summary ──")
print(f"\n {'Item':<38} {'Sc.A':>14} {'Sc.B':>14}")
print(" " + "-" * 68)
print(f" {'Current annual energy cost':<38} £{SCA_COST:>12.2f} £{SCB_OIL_COST:>12.2f}")
print(f" {'Post-Kaluza energy cost':<38} {'Marginal change':>14} £{SCB_ELEC:>12.2f}")
print(f" {'Participation incentive (base)':<38} £{HH_BENEFIT:>11}/yr £{HH_BENEFIT:>11}/yr")
print(f" {'Upfront device cost to household':<38} {'£0':>14} {'£0':>14}")
print(f" {'Net household benefit (base)':<38} £{HH_BENEFIT:>10}/yr+ £{SCB_NET:>12.2f}/yr")
print(f" {'Proposition':<38} {'Positive':>14} {'Weak':>14}")
print(f" {'Adoption driver':<38} {'Incentive':>14} {'Grant+enviro':>14}")
# ── Incentive sensitivity table ───────────────────────────────────────────────
print(f"\n── Sensitivity: Incentive Level ──")
print(f"\n {'Incentive':>12} {'Sc.A net benefit':>18} {'Sc.B net vs oil':>18} {'Sc.B viable?':>12}")
print(" " + "-" * 68)
for inc in [50, 80, 100, 150]:
scb_net = SCB_OIL_COST - SCB_ELEC + inc
print(f" £{inc:>8}/yr £{inc:>14}/yr £{scb_net:>14.2f}/yr "
f"{'✓' if scb_net > 0 else '✗'}")
# ── Household economics before/after comparison ───────────────────────────────
y_max = max(SCA_COST, SCB_ELEC) * 1.35
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
for ax, title, before, after, diff, good in [
(axes[0],
'Sc.A — Existing Storage Heater Household\n(Positive proposition)',
SCA_COST,
SCA_COST - HH_BENEFIT,
HH_BENEFIT,
True),
(axes[1],
'Sc.B — Oil/Coal Household Converting to Storage Heaters\n(Weak proposition at 2017 energy prices)',
SCB_OIL_COST,
SCB_ELEC - HH_BENEFIT,
abs(SCB_NET),
False)
]:
bar_colours = [C5, C3 if good else C2]
ax.bar(['Before\nKaluza', 'After\nKaluza'],
[before, after],
color=bar_colours, alpha=0.85, width=0.4)
ax.text(0, before + y_max * 0.02, f'£{before:.0f}/yr',
ha='center', fontsize=10, fontweight='bold', color='#111111')
ax.text(1, after + y_max * 0.02, f'£{after:.0f}/yr',
ha='center', fontsize=10, fontweight='bold', color='#111111')
diff_label = f'-£{diff:.0f}/yr saving' if good else f'+£{diff:.0f}/yr additional cost'
diff_color = C3 if good else C2
ax.annotate('', xy=(1, after), xytext=(0, before),
arrowprops=dict(arrowstyle='<->', color=diff_color,
lw=1.5, connectionstyle='arc3,rad=0'))
ax.text(0.5, (before + after) / 2, diff_label,
ha='center', fontsize=9, color=diff_color, fontweight='bold',
bbox=dict(boxstyle='round,pad=0.3', facecolor='white',
edgecolor=diff_color, alpha=0.8))
ax.set_ylabel('£ per year')
ax.set_title(title, fontsize=10, fontweight='bold')
ax.set_ylim(0, y_max)
# Footnote on Sc.B after bar
axes[1].text(1, (SCB_ELEC - HH_BENEFIT) - y_max * 0.08,
'£900 E7 cost\n− £80 incentive',
ha='center', fontsize=8, color='#111111', fontweight='bold',
bbox=dict(boxstyle='round,pad=0.2', facecolor='white',
edgecolor='#CCCCCC', alpha=0.9))
fig.suptitle('Figure 9.2: Per-Household Annual Cost — Before and After Kaluza\n'
'(Base case: £80/HH/yr participation incentive, no upfront device cost)',
fontsize=11, fontweight='bold')
plt.tight_layout()
plt.savefig('fig_9_2_household_economics.png', dpi=150, bbox_inches='tight')
plt.show()
AppendixFinancialModel.3 — Household Economics
========================================================================
── Sc.A — Existing Storage Heater Household ──
Annual consumption: 5040 kWh/yr
Ofgem TDCV Aug 2017 (4,200 kWh) × 1.2 Orkney uplift
Night consumption: 2117 kWh/yr (42%)
Day consumption: 2923 kWh/yr (58%)
Night cost: £158.76 (7.5p/kWh E7 off-peak, derived)
Day cost: £467.71 (16.0p/kWh E7 peak, derived)
Current annual cost: £626.47/yr (£52.21/month)
Participation incentive: +£80/yr
Upfront device cost: £0 (Kaluza-funded or grant-covered)
Net household benefit: +£80/yr plus potential scheduling saving
── Sc.B — Oil/Coal Household Converting to Storage Heaters ──
Annual heating demand: 12,000 kWh/yr
Ofgem TDCV Aug 2017 gas medium proxy — flagged in Section 12
Oil price (2017 avg): 39.17p/litre (BEIS QEP Table 4.1.2, 2017)
Calorific value: 10.18 kWh/litre (BEIS conversion factors 2017)
Boiler efficiency: 85% (HHIC industry standard)
Oil effective cost: 4.527p/kWh useful heat (derived)
Current oil cost: £543.21/yr (£45.27/month)
New E7 electricity: £900.00/yr (12,000 kWh × 7.5p)
Participation incentive: +£80/yr
Net vs oil baseline: £-276.79/yr (additional cost)
── Household Economics Summary ──
Item Sc.A Sc.B
--------------------------------------------------------------------
Current annual energy cost £ 626.47 £ 543.21
Post-Kaluza energy cost Marginal change £ 900.00
Participation incentive (base) £ 80/yr £ 80/yr
Upfront device cost to household £0 £0
Net household benefit (base) £ 80/yr+ £ -276.79/yr
Proposition Positive Weak
Adoption driver Incentive Grant+enviro
── Sensitivity: Incentive Level ──
Incentive Sc.A net benefit Sc.B net vs oil Sc.B viable?
--------------------------------------------------------------------
£ 50/yr £ 50/yr £ -306.79/yr ✗
£ 80/yr £ 80/yr £ -276.79/yr ✗
£ 100/yr £ 100/yr £ -256.79/yr ✗
£ 150/yr £ 150/yr £ -206.79/yr ✗
Appendix G: DR Penetration and Commercial Viability¶
Appendix G.1 — DR Penetration¶
In [177]:
"""
AppendixDR.1 — DR Penetration Simulation
==========================================
Answers Case Questions 2 and 3 for Sc.A-only base case:
Q2: How much curtailment can be reduced at different DR penetration levels?
Q3: How many Sc.A households are needed for a given absorption target?
Note: Sc.B is not modelled in the base case and is treated as a future
extension contingent on grant support, tariff innovation, or changes in
relative energy prices (Section 9.3, Section 13).
Sc.A technical ceiling: 32% of total households (3,323 HH).
Commercial base-case upper scenario: 30% penetration.
Inputs:
- data/curtailment_hourly.csv — pre-computed hourly curtailment from
Section 7 notebook, ensuring full methodological consistency.
Method:
For each penetration level:
1. Load fleet-level setpoint-driven curtailment from Section 7 output
2. Determine available DR capacity from enrolled households
3. Apply 70% simultaneous availability factor
4. Apply daily thermal cap constraint per household
5. Compute absorbed MWh as min(available curtailment,
hourly DR capacity,
remaining daily thermal headroom)
Sources:
- Device capacity 2.2 kW/HH: Dimplex Quantum RF input rating
- Daily thermal cap 10 kWh/HH: internal modelling (conservative base case)
Theoretical maximum 15.4 kWh (2.2 kW × 7 hrs Economy 7 window)
Sensitivity at 15 kWh/day in outputs below
- DR availability factor 70%: assumptions register §4, internal modelling
- N_TURBINES 600: assumptions register
- N_HH 10,385: FAQ
"""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
plt.rcParams['figure.dpi'] = 130
plt.rcParams['font.family'] = 'sans-serif'
plt.rcParams['axes.spines.top'] = False
plt.rcParams['axes.spines.right'] = False
SEP = "=" * 72
# ══════════════════════════════════════════════════════════════════════════════
# CONSTANTS
# ══════════════════════════════════════════════════════════════════════════════
N_TURBINES = 600
N_HH = 10_385 # [FAQ]
SCA_SHARE = 0.32 # [OREF EPC 2022]
N_SCA = int(np.floor(N_HH * SCA_SHARE + 0.5)) # 3,323
DEVICE_KW = 2.2 # kW [Dimplex Quantum RF input rating]
DAILY_CAP_KWH = 10.0 # kWh/day — conservative base case [internal modelling]
DAILY_CAP_SENS = 15.0 # kWh/day — sensitivity [theoretical 2.2 kW × 7hr E7 window]
DR_AVAIL_FACTOR = 0.70 # [assumptions register §4]
DR_ADDRESSABLE_MWH = 268_146 # MWh/yr from Section 7
TOTAL_CURTAILED_MWH = 450_431 # MWh/yr from Section 7
# Penetration scenarios — % of TOTAL households [FAQ/coaching/personal notes]
# 30% = commercial base-case upper; 32% = Sc.A technical ceiling
PENETRATIONS = [0.10, 0.15, 0.20, 0.25, 0.30, 0.32]
# ══════════════════════════════════════════════════════════════════════════════
# LOAD DATA
# ══════════════════════════════════════════════════════════════════════════════
print(SEP)
print("AppendixDR.1 — DR Penetration Simulation")
print(SEP)
# Load pre-computed hourly curtailment from Section 7 notebook
# curtailment_hourly.csv uses same methodology as Section 7 fleet summary
df_hourly = pd.read_csv('data/curtailment_hourly.csv',
parse_dates=['Timestamp'],
index_col='Timestamp')
# DR-addressable curtailment = Setpoint_loss_kw
# Hourly data: kW average over one hour is treated as kWh for that hour
df_hourly['dr_addressable_kwh'] = df_hourly['Setpoint_loss_kw']
span_years = (df_hourly.index.max() - df_hourly.index.min()).days / 365.25
print(f"\nLoaded: {len(df_hourly):,} hours over {span_years:.2f} years")
print(f"Date range: {df_hourly.index.min()} to {df_hourly.index.max()}")
computed_dr = df_hourly['dr_addressable_kwh'].sum() / 1000 / span_years
computed_tot = df_hourly['Curtailment_kw'].sum() / 1000 / span_years
print(f"\nFleet curtailment (annualised): {computed_tot:,.0f} MWh/yr")
print(f"DR-addressable (annualised): {computed_dr:,.0f} MWh/yr")
print(f"Section 7 target: {DR_ADDRESSABLE_MWH:,} MWh/yr")
# Consistency check — warn if computed figure deviates materially from Section 7
if abs(computed_dr - DR_ADDRESSABLE_MWH) > 1000:
print(f"WARNING: DR-addressable MWh/yr deviates from Section 7 target "
f"by {abs(computed_dr - DR_ADDRESSABLE_MWH):,.0f} MWh/yr. "
f"Check methodology consistency.")
else:
print(f"Consistency check passed — deviation {abs(computed_dr - DR_ADDRESSABLE_MWH):,.0f} "
f"MWh/yr (<1,000 MWh threshold)")
# ══════════════════════════════════════════════════════════════════════════════
# SIMULATION FUNCTION
# ══════════════════════════════════════════════════════════════════════════════
def simulate_dr(df_hourly, n_enrolled, device_kw, daily_cap_kwh, avail_factor):
"""
Simulate hourly DR absorption for enrolled Sc.A households.
Each hour:
available_hh = n_enrolled × avail_factor
max_hourly_kwh = available_hh × device_kw
max_daily_kwh = available_hh × daily_cap_kwh
absorbed = min(dr_addressable_kwh,
max_hourly_kwh,
remaining_daily_cap)
Daily cap resets at midnight.
"""
available_hh = n_enrolled * avail_factor
max_hourly_kwh = available_hh * device_kw
max_daily_kwh = available_hh * daily_cap_kwh
total_absorbed = 0.0
daily_absorbed = 0.0
current_date = None
for timestamp, row in df_hourly.iterrows():
date = timestamp.date()
if date != current_date:
daily_absorbed = 0.0
current_date = date
available_kwh = row['dr_addressable_kwh']
if available_kwh <= 0:
continue
remaining_cap = max_daily_kwh - daily_absorbed
if remaining_cap <= 0:
continue
absorbed = min(available_kwh, max_hourly_kwh, remaining_cap)
total_absorbed += absorbed
daily_absorbed += absorbed
return total_absorbed
# ══════════════════════════════════════════════════════════════════════════════
# RUN SIMULATION — BASE CASE AND THERMAL CAP SENSITIVITY
# ══════════════════════════════════════════════════════════════════════════════
print(f"\n{SEP}")
print(f"Simulation: device {DEVICE_KW} kW | "
f"availability {DR_AVAIL_FACTOR*100:.0f}% | "
f"base cap {DAILY_CAP_KWH} kWh/day | "
f"sensitivity cap {DAILY_CAP_SENS} kWh/day")
print(SEP)
results = []
for pen in PENETRATIONS:
n_enrolled = min(int(np.floor(N_HH * pen + 0.5)), N_SCA)
absorbed_base = simulate_dr(df_hourly, n_enrolled,
DEVICE_KW, DAILY_CAP_KWH, DR_AVAIL_FACTOR)
absorbed_sens = simulate_dr(df_hourly, n_enrolled,
DEVICE_KW, DAILY_CAP_SENS, DR_AVAIL_FACTOR)
mwh_base = (absorbed_base / 1000) / span_years
mwh_sens = (absorbed_sens / 1000) / span_years
flex_mw = n_enrolled * DR_AVAIL_FACTOR * DEVICE_KW / 1000
pct_addr_base = 100 * mwh_base / DR_ADDRESSABLE_MWH
pct_tot_base = 100 * mwh_base / TOTAL_CURTAILED_MWH
is_ceiling = (n_enrolled == N_SCA)
pen_label = f"{pen*100:.0f}%" + (" (Sc.A ceiling)" if is_ceiling else "")
results.append({
'Penetration (%)': pen_label,
'HH enrolled': n_enrolled,
'Flex capacity (MW)': round(flex_mw, 2),
'Absorbed — base (MWh/yr)': round(mwh_base),
'Absorbed — sens (MWh/yr)': round(mwh_sens),
'% of DR-addressable (base)': round(pct_addr_base, 2),
'% of total curtailment': round(pct_tot_base, 2),
})
print(f" {pen_label:<20}: "
f"{n_enrolled:,} HH | {flex_mw:.2f} MW | "
f"base {mwh_base:,.0f} MWh/yr ({pct_addr_base:.2f}% addr) | "
f"sens {mwh_sens:,.0f} MWh/yr")
print("\n * = Sc.A technical ceiling (32%, 3,323 HH)")
df_results = pd.DataFrame(results)
print(f"\n{SEP}")
print("PENETRATION SIMULATION RESULTS")
print(SEP)
print(df_results.to_string(index=False))
df_results.to_csv('dr_penetration_results.csv', index=False)
# ══════════════════════════════════════════════════════════════════════════════
# CASE QUESTION 3 — HOUSEHOLDS REQUIRED BY ABSORPTION TARGET
# ══════════════════════════════════════════════════════════════════════════════
print(f"\n{SEP}")
print("AppendixDR.1 — Case Question 3: Households Required by Absorption Target")
print(SEP)
# Base case absorbed values including (0,0) origin for correct interpolation
absorbed_vals = [0] + [r['Absorbed — base (MWh/yr)'] for r in results]
enrolled_vals = [0] + [r['HH enrolled'] for r in results]
penetration_vals = [0.0] + [float(r['Penetration (%)'].split('%')[0])/100
for r in results]
target_reductions_pct = [0.1, 0.25, 0.5, 1.0, 2.0]
print(f"\n {'Target':>22} {'MWh/yr':>10} {'HH needed':>10} {'Penetration':>12}")
print(" " + "-" * 60)
for target_pct in target_reductions_pct:
target_mwh = DR_ADDRESSABLE_MWH * target_pct / 100
if target_mwh <= 0:
continue
elif target_mwh >= absorbed_vals[-1]:
print(f" {target_pct:.1f}% of addressable {target_mwh:>10,.0f} "
f"{enrolled_vals[-1]:>10,} >Sc.A ceiling")
continue
else:
for i in range(len(absorbed_vals)-1):
if absorbed_vals[i] <= target_mwh <= absorbed_vals[i+1]:
frac = ((target_mwh - absorbed_vals[i]) /
(absorbed_vals[i+1] - absorbed_vals[i]))
hh = int(enrolled_vals[i] + frac*(enrolled_vals[i+1]-enrolled_vals[i]))
pen = penetration_vals[i] + frac*(penetration_vals[i+1]-penetration_vals[i])
print(f" {target_pct:.1f}% of addressable {target_mwh:>10,.0f} "
f"{hh:>10,} ~{pen*100:.1f}%")
break
# ══════════════════════════════════════════════════════════════════════════════
# VISUALISATIONS
# ══════════════════════════════════════════════════════════════════════════════
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
pct_labels = df_results['Penetration (%)'].tolist()
abs_base = df_results['Absorbed — base (MWh/yr)'].tolist()
pct_addr = df_results['% of DR-addressable (base)'].tolist()
# ─────────────────────────────────────────────────────────────────────────────
# Figure DR.1a — Curtailment absorbed
# ─────────────────────────────────────────────────────────────────────────────
axes[0].bar(pct_labels, abs_base, color=C1, alpha=0.85)
axes[0].set_xlabel('DR Penetration (% of Total Households)')
axes[0].set_ylabel('Curtailment Absorbed (MWh/year)')
axes[0].set_title('Figure G.1a: Curtailment Absorbed by Penetration Level')
for i, v in enumerate(abs_base):
axes[0].text(
i,
v + max(abs_base) * 0.015,
f'{v:,.0f}',
ha='center',
fontsize=8,
color='#111111'
)
# ─────────────────────────────────────────────────────────────────────────────
# Figure DR.1b — % of DR-addressable pool
# ─────────────────────────────────────────────────────────────────────────────
axes[1].bar(pct_labels, pct_addr, color=C3, alpha=0.85)
axes[1].set_xlabel('DR Penetration (% of Total Households)')
axes[1].set_ylabel('Share of DR-Addressable Curtailment (%)')
axes[1].set_title('Figure G.1b: Absorption as a Share of the\nDR-Addressable Curtailment Pool')
for i, v in enumerate(pct_addr):
axes[1].text(
i,
v + max(pct_addr) * 0.03,
f'{v:.2f}%',
ha='center',
fontsize=8,
color='#111111'
)
plt.tight_layout()
plt.savefig('fig_dr_penetration.png', dpi=150, bbox_inches='tight')
plt.show()
# ══════════════════════════════════════════════════════════════════════════════
# KEY OUTPUTS SUMMARY
# ══════════════════════════════════════════════════════════════════════════════
print(f"\n{SEP}")
print("KEY OUTPUTS — AppendixDR.1")
print(SEP)
max_res = next(r for r in results if r['Penetration (%)'] == '30%')
ceil_res = next(r for r in results if 'ceiling' in r['Penetration (%)'])
print(f"\nDR-addressable pool (Section 7): {DR_ADDRESSABLE_MWH:,} MWh/yr")
print(f"Total curtailment (Section 7): {TOTAL_CURTAILED_MWH:,} MWh/yr")
print(f"\nBase-case upper scenario (30%):")
print(f" HH enrolled: {max_res['HH enrolled']:,}")
print(f" Flexible demand capacity: {max_res['Flex capacity (MW)']:.2f} MW")
print(f" Absorbed (base): {max_res['Absorbed — base (MWh/yr)']:,} MWh/yr")
print(f" Absorbed (sensitivity): {max_res['Absorbed — sens (MWh/yr)']:,} MWh/yr")
print(f" % of DR-addressable (base): {max_res['% of DR-addressable (base)']:.2f}%")
print(f"\nSc.A technical ceiling (32%, {N_SCA:,} HH):")
print(f" Absorbed (base): {ceil_res['Absorbed — base (MWh/yr)']:,} MWh/yr")
print(f" Absorbed (sensitivity): {ceil_res['Absorbed — sens (MWh/yr)']:,} MWh/yr")
========================================================================
AppendixDR.1 — DR Penetration Simulation
========================================================================
Loaded: 17,852 hours over 2.63 years
Date range: 2015-05-28 00:00:00 to 2018-01-11 06:00:00
Fleet curtailment (annualised): 450,857 MWh/yr
DR-addressable (annualised): 268,400 MWh/yr
Section 7 target: 268,146 MWh/yr
Consistency check passed — deviation 254 MWh/yr (<1,000 MWh threshold)
========================================================================
Simulation: device 2.2 kW | availability 70% | base cap 10.0 kWh/day | sensitivity cap 15.0 kWh/day
========================================================================
10% : 1,039 HH | 1.60 MW | base 646 MWh/yr (0.24% addr) | sens 880 MWh/yr
15% : 1,558 HH | 2.40 MW | base 963 MWh/yr (0.36% addr) | sens 1,313 MWh/yr
20% : 2,077 HH | 3.20 MW | base 1,277 MWh/yr (0.48% addr) | sens 1,743 MWh/yr
25% : 2,596 HH | 4.00 MW | base 1,589 MWh/yr (0.59% addr) | sens 2,169 MWh/yr
30% : 3,116 HH | 4.80 MW | base 1,900 MWh/yr (0.71% addr) | sens 2,594 MWh/yr
32% (Sc.A ceiling) : 3,323 HH | 5.12 MW | base 2,023 MWh/yr (0.75% addr) | sens 2,762 MWh/yr
* = Sc.A technical ceiling (32%, 3,323 HH)
========================================================================
PENETRATION SIMULATION RESULTS
========================================================================
Penetration (%) HH enrolled Flex capacity (MW) Absorbed — base (MWh/yr) Absorbed — sens (MWh/yr) % of DR-addressable (base) % of total curtailment
10% 1039 1.60 646 880 0.24 0.14
15% 1558 2.40 963 1313 0.36 0.21
20% 2077 3.20 1277 1743 0.48 0.28
25% 2596 4.00 1589 2169 0.59 0.35
30% 3116 4.80 1900 2594 0.71 0.42
32% (Sc.A ceiling) 3323 5.12 2023 2762 0.75 0.45
========================================================================
AppendixDR.1 — Case Question 3: Households Required by Absorption Target
========================================================================
Target MWh/yr HH needed Penetration
------------------------------------------------------------
0.1% of addressable 268 431 ~4.2%
0.2% of addressable 670 1,078 ~10.4%
0.5% of addressable 1,341 2,183 ~21.0%
1.0% of addressable 2,681 3,323 >Sc.A ceiling
2.0% of addressable 5,363 3,323 >Sc.A ceiling
======================================================================== KEY OUTPUTS — AppendixDR.1 ======================================================================== DR-addressable pool (Section 7): 268,146 MWh/yr Total curtailment (Section 7): 450,431 MWh/yr Base-case upper scenario (30%): HH enrolled: 3,116 Flexible demand capacity: 4.80 MW Absorbed (base): 1,900 MWh/yr Absorbed (sensitivity): 2,594 MWh/yr % of DR-addressable (base): 0.71% Sc.A technical ceiling (32%, 3,323 HH): Absorbed (base): 2,023 MWh/yr Absorbed (sensitivity): 2,762 MWh/yr
Appendix G.2 - Commercial Viability Analysis¶
In [164]:
# ============================================================
# AppendixDR.2 — Commercial Viability Analysis
# Corresponds to Section 10.2 of the main report.
# Defines kaluza_pnl() and reads dr_penetration_results.csv
# from AppendixDR.1.
# ============================================================
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
SEP = "=" * 72
dr = pd.read_csv('dr_penetration_results.csv')
PENETRATIONS_DR = dr['Penetration (%)'].tolist()
ABSORBED_MWH = dr['Absorbed — base (MWh/yr)'].tolist()
HH_ENROLLED = dr['HH enrolled'].tolist()
print(SEP)
print("AppendixDR.2 — Commercial Viability Analysis")
print("Sc.A only | Accounting P&L and cash-flow view")
print(SEP)
def eligible_split(n_enrolled):
n_elig = round(n_enrolled * FUEL_POV)
n_nonelig = n_enrolled - n_elig
return n_elig, n_nonelig
def kaluza_pnl(n_enrolled, absorbed_mwh, ssen_rate, gen_fee,
extra_subsidy=0.0, incentive=None):
"""
Annual steady-state accounting P&L — Sc.A only.
Capex and onboarding are annualised over DEVICE_LIFE.
"""
if incentive is None:
incentive = KALUZA_INC
ssen_rev = absorbed_mwh * ssen_rate
gen_rev = absorbed_mwh * gen_fee
total_rev = ssen_rev + gen_rev
_, n_nonelig = eligible_split(n_enrolled)
capex_net = n_nonelig * CAPEX_A * (1 - extra_subsidy)
acq_net = n_enrolled * ONE_OFF
capex_yr = (capex_net + acq_net) / DEVICE_LIFE
opex = n_enrolled * (MAINT + incentive) + FIXED_OPEX
total_cost = capex_yr + opex
profit = total_rev - total_cost
return {
'n_enrolled': n_enrolled,
'ssen_rev': round(ssen_rev),
'gen_rev': round(gen_rev),
'total_rev': round(total_rev),
'capex_yr': round(capex_yr),
'opex': round(opex),
'total_cost': round(total_cost),
'profit': round(profit),
'viable': profit > 0,
}
print(f"\n── Annual Accounting P&L by Penetration ──")
print(f"SSEN £{SSEN_BASE}/MWh | Gen fee £{GEN_FEE_BASE}/MWh | 0% extra subsidy\n")
print(f"{'Penetration':>20} {'HH':>6} {'MWh/yr':>8} {'Revenue':>12} "
f"{'Cost':>12} {'Profit':>12} {'Viable':>7}")
print("-" * 80)
for pen, mwh, hh in zip(PENETRATIONS_DR, ABSORBED_MWH, HH_ENROLLED):
r = kaluza_pnl(hh, mwh, SSEN_BASE, GEN_FEE_BASE, 0.0)
print(f" {pen:>18} {hh:>6,} {mwh:>8,} "
f"£{r['total_rev']:>10,} £{r['total_cost']:>10,} "
f"£{r['profit']:>10,} {'✓' if r['viable'] else '✗':>6}")
print(f"\n── Break-even Penetration — Accounting P&L ──\n")
found = False
for pen, mwh, hh in zip(PENETRATIONS_DR, ABSORBED_MWH, HH_ENROLLED):
if 'ceiling' in str(pen):
continue
r = kaluza_pnl(hh, mwh, SSEN_BASE, GEN_FEE_BASE, 0.0)
if r['viable'] and not found:
print(f" First break-even at {pen} penetration ({hh:,} households)")
print(f" Annual accounting profit: £{r['profit']:,}")
print(f" Absorbed MWh: {mwh:,} MWh/yr")
found = True
if not found:
print(" Does not break even within Sc.A pool under base-case assumptions")
print(f"\n── Break-even Penetration by Subsidy Scenario — Accounting P&L ──\n")
for sub in [0.0, 0.25, 0.50, 0.75]:
found = False
for pen, mwh, hh in zip(PENETRATIONS_DR, ABSORBED_MWH, HH_ENROLLED):
if 'ceiling' in str(pen):
continue
r = kaluza_pnl(hh, mwh, SSEN_BASE, GEN_FEE_BASE, sub)
if r['viable'] and not found:
print(f" {sub*100:.0f}% subsidy: first break-even at {pen} "
f"({hh:,} HH | annual profit £{r['profit']:,})")
found = True
if not found:
print(f" {sub*100:.0f}% subsidy: does not break even within Sc.A pool")
print(f"\n── Minimum SSEN Rate for Accounting Profitability ──\n")
for pen, mwh, hh in zip(PENETRATIONS_DR, ABSORBED_MWH, HH_ENROLLED):
if 'ceiling' in str(pen):
continue
found_rate = False
for s in range(50, 600, 5):
r = kaluza_pnl(hh, mwh, s, GEN_FEE_BASE, 0.0)
if r['viable']:
gap = abs(SSEN_BASE - s)
pos = 'below' if s < SSEN_BASE else 'at' if s == SSEN_BASE else 'above'
print(f" {pen}: needs SSEN ≥ £{s}/MWh "
f"£{gap} {pos} base-case £{SSEN_BASE}/MWh")
found_rate = True
break
if not found_rate:
print(f" {pen}: no break-even SSEN rate found within tested range")
print(f"\n{SEP}")
print("5-Year Cash-Flow Projection — Base Case")
print(SEP)
RAMP = {2018: 0.30, 2019: 0.70, 2020: 1.00, 2021: 1.00, 2022: 1.00}
ADDITIONS = {2018: 0.30, 2019: 0.40, 2020: 0.30, 2021: 0.00, 2022: 0.00}
DISC = 0.10
print(f"\n── 5-Year Summary Across Penetration Targets ──\n")
print(f"{'Target':>20} {'HH':>6} {'Cum Cash Flow':>16} "
f"{'NPV Cash Flow':>16} {'Positive Yr5':>14}")
print("-" * 82)
proj_30 = None
for pen, mwh, hh in zip(PENETRATIONS_DR, ABSORBED_MWH, HH_ENROLLED):
if 'ceiling' in str(pen):
continue
_, n_nonelig = eligible_split(hh)
total_initial_capex = n_nonelig * CAPEX_A + hh * ONE_OFF
cumulative_cf = 0.0
npv_cf = 0.0
rows = []
for i, yr in enumerate([2018, 2019, 2020, 2021, 2022], start=1):
ramp = RAMP[yr]
add = ADDITIONS[yr]
revenue = mwh * ramp * (SSEN_BASE + GEN_FEE_BASE)
opex = (hh * OPEX_PER_HH + FIXED_OPEX) * ramp
rollout_capex = total_initial_capex * add
cash_flow = revenue - opex - rollout_capex
cumulative_cf += cash_flow
discounted_cf = cash_flow / (1 + DISC) ** i
npv_cf += discounted_cf
rows.append({
'Year': yr,
'Ramp': f"{ramp*100:.0f}%",
'HH active': int(hh * ramp),
'MWh': round(mwh * ramp),
'Revenue (£)': round(revenue),
'Opex (£)': round(opex),
'Rollout capex (£)': round(rollout_capex),
'Cash flow (£)': round(cash_flow),
'Cumulative cash flow (£)': round(cumulative_cf),
'Discounted cash flow (£)': round(discounted_cf),
})
print(f" {pen:>18} {hh:>6,} £{round(cumulative_cf):>14,} "
f"£{round(npv_cf):>14,} {'YES' if cumulative_cf > 0 else 'NO':>14}")
if pen == '30%':
proj_30 = (pd.DataFrame(rows), round(npv_cf), hh)
if proj_30:
df_proj, npv_30, hh_30 = proj_30
print(f"\nDetailed 5-year cash-flow projection — 30% penetration ({hh_30:,} HH):\n")
print(df_proj.to_string(index=False))
print(f"\n Cumulative cash flow: £{df_proj['Cumulative cash flow (£)'].iloc[-1]:,}")
print(f" NPV cash flow (10%): £{npv_30:,}")
# ── Figure 10.1: 5-Year Cumulative Cash Flow ─────────────────────────────────
RAMP = {2018: 0.30, 2019: 0.70, 2020: 1.00, 2021: 1.00, 2022: 1.00}
ADDITIONS = {2018: 0.30, 2019: 0.40, 2020: 0.30, 2021: 0.00, 2022: 0.00}
DISC = 0.10
YEARS = [2018, 2019, 2020, 2021, 2022]
YR_LABELS = [str(y) for y in YEARS]
# HSO palette for penetration lines
colors_pen = {
'10%': C2,
'15%': C4,
'20%': C3,
'25%': C5,
'30%': C1,
}
plot_items = [
(p, m, h) for p, m, h in zip(PENETRATIONS_DR, ABSORBED_MWH, HH_ENROLLED)
if 'ceiling' not in str(p)
]
fig, ax = plt.subplots(figsize=(10, 5))
for pen, mwh, hh in plot_items:
col = colors_pen.get(pen, C5)
_, n_nonelig = eligible_split(hh)
total_capex = n_nonelig * CAPEX_A + hh * ONE_OFF
cum_cf = 0.0
cum_vals = []
for yr in YEARS:
ramp = RAMP[yr]
add = ADDITIONS[yr]
rev = mwh * ramp * (SSEN_BASE + GEN_FEE_BASE)
opex = (hh * OPEX_PER_HH + FIXED_OPEX) * ramp
cap = total_capex * add
cum_cf += rev - opex - cap
cum_vals.append(cum_cf / 1000)
ax.plot(YR_LABELS, cum_vals, marker='o', label=pen,
color=col, linewidth=2)
ax.axhline(0, color='#444444', linewidth=1, linestyle='--', alpha=0.7)
ax.set_xlabel('Year')
ax.set_ylabel('Cumulative cash flow (£k)')
ax.set_title(
'Figure 10.2: Five-Year Cumulative Cash Flow by Sc.A Penetration Scenario',
fontsize=12,
fontweight='bold',
pad=20
)
ax.text(
0.5, 1.02,
'Base-case assumptions: SSEN £250/MWh, wind farmer fee £25/MWh, 0% additional subsidy',
transform=ax.transAxes,
ha='center',
fontsize=10,
color='#555555'
)
ax.legend(title='Penetration', fontsize=8)
ax.yaxis.set_major_formatter(mticker.FuncFormatter(lambda x, _: f'£{x:,.0f}k'))
plt.tight_layout()
plt.savefig('fig_10_1_cumulative_cashflow.png', dpi=150, bbox_inches='tight')
plt.show()
# ── Figure 10.2: Revenue vs Cost vs Profit by Penetration ────────────────────
fig, ax = plt.subplots(figsize=(12, 5))
x = np.arange(len(plot_items))
w = 0.25
revs = []
costs = []
profits = []
labels = []
for pen, mwh, hh in plot_items:
r = kaluza_pnl(hh, mwh, SSEN_BASE, GEN_FEE_BASE, 0.0)
revs.append(r['total_rev'] / 1000)
costs.append(r['total_cost'] / 1000)
profits.append(r['profit'] / 1000)
labels.append(pen)
ax.bar(x - w, revs, w, label='Revenue', color=C3, alpha=0.85)
ax.bar(x, costs, w, label='Total cost', color=C2, alpha=0.85)
ax.bar(x + w, profits, w, label='Profit/Loss',color=C1, alpha=0.85)
# Annotate profit/loss bars
for i, v in enumerate(profits):
ax.text(i + w, v + (max(revs) * 0.01 if v >= 0 else -max(revs) * 0.03),
f'£{v:,.0f}k', ha='center', fontsize=7,
color=C1, fontweight='bold')
ax.axhline(0, color='#444444', linewidth=0.8)
ax.set_xticks(x)
ax.set_xticklabels(labels)
ax.set_xlabel('Sc.A Penetration (% of total households)')
ax.set_ylabel('£k per year (steady state)')
ax.set_title(
'Figure 10.1: Annual Revenue, Cost and Accounting Profit by Sc.A Penetration Scenario',
fontsize=12,
fontweight='bold',
pad=20
)
ax.text(
0.5, 1.02,
f'Steady-state operation under base-case assumptions '
f'(SSEN £{SSEN_BASE}/MWh, wind farmer fee £{GEN_FEE_BASE}/MWh)',
transform=ax.transAxes,
ha='center',
fontsize=10,
color='#555555'
)
ax.legend(fontsize=9)
ax.yaxis.set_major_formatter(mticker.FuncFormatter(lambda x, _: f'£{x:,.0f}k'))
plt.tight_layout()
plt.savefig('fig_10_2_revenue_cost_profit.png', dpi=150, bbox_inches='tight')
plt.show()
========================================================================
AppendixDR.2 — Commercial Viability Analysis
Sc.A only | Accounting P&L and cash-flow view
========================================================================
── Annual Accounting P&L by Penetration ──
SSEN £250/MWh | Gen fee £25/MWh | 0% extra subsidy
Penetration HH MWh/yr Revenue Cost Profit Viable
--------------------------------------------------------------------------------
10% 1,039 646 £ 177,650 £ 180,551 £ -2,901 ✗
15% 1,558 963 £ 264,825 £ 228,234 £ 36,590 ✓
20% 2,077 1,277 £ 351,175 £ 275,918 £ 75,257 ✓
25% 2,596 1,589 £ 436,975 £ 323,602 £ 113,374 ✓
30% 3,116 1,900 £ 522,500 £ 371,388 £ 151,112 ✓
32% (Sc.A ceiling) 3,323 2,023 £ 556,325 £ 390,407 £ 165,918 ✓
── Break-even Penetration — Accounting P&L ──
First break-even at 15% penetration (1,558 households)
Annual accounting profit: £36,590
Absorbed MWh: 963 MWh/yr
── Break-even Penetration by Subsidy Scenario — Accounting P&L ──
0% subsidy: first break-even at 15% (1,558 HH | annual profit £36,590)
25% subsidy: first break-even at 15% (1,558 HH | annual profit £39,656)
50% subsidy: first break-even at 10% (1,039 HH | annual profit £1,186)
75% subsidy: first break-even at 10% (1,039 HH | annual profit £3,230)
── Minimum SSEN Rate for Accounting Profitability ──
10%: needs SSEN ≥ £255/MWh £5 above base-case £250/MWh
15%: needs SSEN ≥ £215/MWh £35 below base-case £250/MWh
20%: needs SSEN ≥ £195/MWh £55 below base-case £250/MWh
25%: needs SSEN ≥ £180/MWh £70 below base-case £250/MWh
30%: needs SSEN ≥ £175/MWh £75 below base-case £250/MWh
========================================================================
5-Year Cash-Flow Projection — Base Case
========================================================================
── 5-Year Summary Across Penetration Targets ──
Target HH Cum Cash Flow NPV Cash Flow Positive Yr5
----------------------------------------------------------------------------------
10% 1,039 £ -118,564 £ -102,529 NO
15% 1,558 £ -14,058 £ -34,778 NO
20% 2,077 £ 87,148 £ 30,575 YES
25% 2,596 £ 186,154 £ 94,329 YES
30% 3,116 £ 283,459 £ 156,819 YES
Detailed 5-year cash-flow projection — 30% penetration (3,116 HH):
Year Ramp HH active MWh Revenue (£) Opex (£) Rollout capex (£) Cash flow (£) Cumulative cash flow (£) Discounted cash flow (£)
2018 30% 934 570 156750 99379 144446 -87075 -87075 -79159
2019 70% 2181 1330 365750 231885 192594 -58729 -145804 -48536
2020 100% 3116 1900 522500 331264 144446 46790 -99013 35154
2021 100% 3116 1900 522500 331264 0 191236 92223 130617
2022 100% 3116 1900 522500 331264 0 191236 283459 118743
Cumulative cash flow: £283,459
NPV cash flow (10%): £156,819
Appendix G.3 - Sensitivity Analysis¶
In [ ]:
# ============================================================
# AppendixDR.3 — Sensitivity Analysis
# Corresponds to Section 10.3 of the main report.
# Assumes all constants from AppendixFinancialModel.1 and .2
# and kaluza_pnl() from AppendixDR.2 already defined.
# Reads dr_penetration_results.csv from AppendixDR.1.
# ============================================================
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
SEP = "=" * 72
# ── Load DR simulation results ────────────────────────────────────────────────
dr = pd.read_csv('dr_penetration_results.csv')
PENETRATIONS_DR = dr['Penetration (%)'].tolist()
ABSORBED_MWH = dr['Absorbed — base (MWh/yr)'].tolist()
HH_ENROLLED = dr['HH enrolled'].tolist()
# 30% base case figures
idx_30 = next(i for i, p in enumerate(PENETRATIONS_DR) if p == '30%')
HH_30 = HH_ENROLLED[idx_30]
MWH_30 = ABSORBED_MWH[idx_30]
MWH_30_SENS = dr[dr['Penetration (%)'] == '30%']['Absorbed — sens (MWh/yr)'].values[0]
pct_labels = [p for p in PENETRATIONS_DR if 'ceiling' not in str(p)]
mwh_vals = [m for p, m in zip(PENETRATIONS_DR, ABSORBED_MWH) if 'ceiling' not in str(p)]
hh_vals = [h for p, h in zip(PENETRATIONS_DR, HH_ENROLLED) if 'ceiling' not in str(p)]
print(SEP)
print("AppendixDR.3 — Sensitivity Analysis")
print("30% Sc.A penetration base case unless stated")
print(SEP)
base_r = kaluza_pnl(HH_30, MWH_30, SSEN_BASE, GEN_FEE_BASE, 0.0)
base_pnl = base_r['profit']
print(f"\nBase case annual profit (30%, SSEN £{SSEN_BASE}/MWh, "
f"wind farmer fee £{GEN_FEE_BASE}/MWh, 0% subsidy): £{base_pnl:,}")
# ══════════════════════════════════════════════════════════════════════════════
# SSEN PRICE SENSITIVITY
# ══════════════════════════════════════════════════════════════════════════════
print(f"\n── SSEN Price Sensitivity (wind farmer fee £{GEN_FEE_BASE}/MWh, 0% subsidy) ──\n")
print(f"{'Penetration':>20}", end="")
for s in SSEN_SCENS:
print(f" £{s}/MWh", end="")
print()
print("-" * 80)
for pen, mwh, hh in zip(pct_labels, mwh_vals, hh_vals):
print(f" {pen:>18}", end="")
for s in SSEN_SCENS:
r = kaluza_pnl(hh, mwh, s, GEN_FEE_BASE, 0.0)
m = '✓' if r['viable'] else '✗'
print(f" £{r['profit']:>8,}{m}", end="")
print()
# ══════════════════════════════════════════════════════════════════════════════
# WIND FARMER FEE SENSITIVITY
# ══════════════════════════════════════════════════════════════════════════════
print(f"\n── Wind Farmer Fee Sensitivity (SSEN £{SSEN_BASE}/MWh, 0% subsidy) ──\n")
print(f"{'Penetration':>20}", end="")
for gf in GEN_FEE_SCENS:
print(f" £{gf}/MWh fee", end="")
print()
print("-" * 70)
for pen, mwh, hh in zip(pct_labels, mwh_vals, hh_vals):
print(f" {pen:>18}", end="")
for gf in GEN_FEE_SCENS:
r = kaluza_pnl(hh, mwh, SSEN_BASE, gf, 0.0)
m = '✓' if r['viable'] else '✗'
print(f" £{r['profit']:>8,}{m}", end="")
print()
# ══════════════════════════════════════════════════════════════════════════════
# GOVERNMENT SUBSIDY SCENARIOS
# ══════════════════════════════════════════════════════════════════════════════
print(f"\n── Government Subsidy Scenarios (SSEN £{SSEN_BASE}/MWh, wind farmer fee £{GEN_FEE_BASE}/MWh) ──\n")
print(f"{'Penetration':>20} {'0% subsidy':>14} {'25% subsidy':>14} "
f"{'50% subsidy':>14} {'75% subsidy':>14}")
print("-" * 82)
for pen, mwh, hh in zip(pct_labels, mwh_vals, hh_vals):
profits = [kaluza_pnl(hh, mwh, SSEN_BASE, GEN_FEE_BASE, sub)['profit']
for sub in [0.0, 0.25, 0.50, 0.75]]
markers = ['✓' if p > 0 else '✗' for p in profits]
print(f" {pen:>18} "
f" £{profits[0]:>10,}{markers[0]} "
f" £{profits[1]:>10,}{markers[1]} "
f" £{profits[2]:>10,}{markers[2]} "
f" £{profits[3]:>10,}{markers[3]}")
# ══════════════════════════════════════════════════════════════════════════════
# INCENTIVE SENSITIVITY
# ══════════════════════════════════════════════════════════════════════════════
KALUZA_INC_HI = round(HH_BENEFIT_HI * 0.60) # £60/HH — upper bound Kaluza portion
print(f"\n── Participation Incentive Sensitivity (30%, SSEN £{SSEN_BASE}/MWh) ──\n")
for inc, label in [(KALUZA_INC, 'base £48/HH'), (KALUZA_INC_HI, 'upper £60/HH')]:
r = kaluza_pnl(HH_30, MWH_30, SSEN_BASE, GEN_FEE_BASE, 0.0, inc)
print(f" £{inc}/HH/yr Kaluza incentive ({label}): P&L £{r['profit']:,} "
f"{'✓' if r['viable'] else '✗'}")
# ══════════════════════════════════════════════════════════════════════════════
# THERMAL CAP SENSITIVITY
# ══════════════════════════════════════════════════════════════════════════════
print(f"\n── Thermal Cap Sensitivity (30%, SSEN £{SSEN_BASE}/MWh) ──\n")
for mwh_val, label in [(MWH_30, f'Base (10 kWh/day): {MWH_30:,} MWh/yr'),
(MWH_30_SENS, f'Sensitivity (15 kWh/day): {MWH_30_SENS:,} MWh/yr')]:
r = kaluza_pnl(HH_30, mwh_val, SSEN_BASE, GEN_FEE_BASE, 0.0)
print(f" {label} → P&L £{r['profit']:,} {'✓' if r['viable'] else '✗'}")
# ══════════════════════════════════════════════════════════════════════════════
# TORNADO CHART
# ══════════════════════════════════════════════════════════════════════════════
tornado = {
f'SSEN £{min(SSEN_SCENS)}/MWh (downside)': kaluza_pnl(HH_30, MWH_30, min(SSEN_SCENS), GEN_FEE_BASE, 0.0)['profit'] - base_pnl,
f'SSEN £{max(SSEN_SCENS)}/MWh (upside)': kaluza_pnl(HH_30, MWH_30, max(SSEN_SCENS), GEN_FEE_BASE, 0.0)['profit'] - base_pnl,
f'Wind farmer fee £{min(GEN_FEE_SCENS)}/MWh': kaluza_pnl(HH_30, MWH_30, SSEN_BASE, min(GEN_FEE_SCENS), 0.0)['profit'] - base_pnl,
f'Wind farmer fee £{max(GEN_FEE_SCENS)}/MWh': kaluza_pnl(HH_30, MWH_30, SSEN_BASE, max(GEN_FEE_SCENS), 0.0)['profit'] - base_pnl,
f'Government subsidy 25%': kaluza_pnl(HH_30, MWH_30, SSEN_BASE, GEN_FEE_BASE, 0.25)['profit'] - base_pnl,
f'Government subsidy 75%': kaluza_pnl(HH_30, MWH_30, SSEN_BASE, GEN_FEE_BASE, 0.75)['profit'] - base_pnl,
f'Incentive £{KALUZA_INC_HI}/HH Kaluza (upper)': kaluza_pnl(HH_30, MWH_30, SSEN_BASE, GEN_FEE_BASE, 0.0, KALUZA_INC_HI)['profit'] - base_pnl,
}
sorted_t = sorted(tornado.items(), key=lambda x: abs(x[1]), reverse=True)
print(f"\n── Tornado — 30% Penetration | Base profit: £{base_pnl:,} ──\n")
for label, delta in sorted_t:
d = '▲' if delta >= 0 else '▼'
print(f" {d} {label:<46} £{delta:+,}")
# ── Figure 10.3: Sensitivity Tornado ─────────────────────────────────────────
fig, ax = plt.subplots(figsize=(10, 6))
t_labels = [x[0] for x in sorted_t]
t_vals = [x[1]/1000 for x in sorted_t]
t_colors = [C3 if v >= 0 else C2 for v in t_vals]
ax.barh(t_labels, t_vals, color=t_colors, alpha=0.85)
ax.axvline(0, color='#444444', linewidth=0.8)
ax.set_title(
'Figure 10.3: Sensitivity of Annual Profit to Key Assumptions',
fontsize=16,
fontweight='bold',
pad=28
)
ax.text(
0.5, 1.02,
'30% Sc.A penetration (SSEN £250/MWh, wind farmer fee £25/MWh, 0% additional subsidy)',
transform=ax.transAxes,
ha='center',
va='bottom',
fontsize=11,
color='dimgray'
)
ax.set_xlabel('Impact on annual profit relative to base case (£k)')
ax.set_ylabel('Sensitivity variable')
ax.xaxis.set_major_formatter(mticker.FuncFormatter(lambda x, _: f'£{x:,.0f}k'))
plt.tight_layout()
plt.savefig('fig_10_3_tornado.png', dpi=150, bbox_inches='tight')
plt.show()
========================================================================
AppendixDR.3 — Sensitivity Analysis
30% Sc.A penetration base case unless stated
========================================================================
Base case annual profit (30%, SSEN £250/MWh, wind farmer fee £25/MWh, 0% subsidy): £151,112
── SSEN Price Sensitivity (wind farmer fee £25/MWh, 0% subsidy) ──
Penetration £150/MWh £200/MWh £250/MWh £300/MWh £400/MWh
--------------------------------------------------------------------------------
10% £ -67,501✗ £ -35,201✗ £ -2,901✗ £ 29,399✓ £ 93,999✓
15% £ -59,710✗ £ -11,560✗ £ 36,590✓ £ 84,740✓ £ 181,040✓
20% £ -52,443✗ £ 11,407✓ £ 75,257✓ £ 139,107✓ £ 266,807✓
25% £ -45,526✗ £ 33,924✓ £ 113,374✓ £ 192,824✓ £ 351,724✓
30% £ -38,888✗ £ 56,112✓ £ 151,112✓ £ 246,112✓ £ 436,112✓
── Wind Farmer Fee Sensitivity (SSEN £250/MWh, 0% subsidy) ──
Penetration £0/MWh fee £15/MWh fee £25/MWh fee £40/MWh fee
----------------------------------------------------------------------
10% £ -19,051✗ £ -9,361✗ £ -2,901✗ £ 6,789✓
15% £ 12,516✓ £ 26,960✓ £ 36,590✓ £ 51,036✓
20% £ 43,332✓ £ 62,487✓ £ 75,257✓ £ 94,412✓
25% £ 73,648✓ £ 97,484✓ £ 113,374✓ £ 137,208✓
30% £ 103,612✓ £ 132,112✓ £ 151,112✓ £ 179,612✓
── Government Subsidy Scenarios (SSEN £250/MWh, wind farmer fee £25/MWh) ──
Penetration 0% subsidy 25% subsidy 50% subsidy 75% subsidy
----------------------------------------------------------------------------------
10% £ -2,901✗ £ -857✗ £ 1,186✓ £ 3,230✓
15% £ 36,590✓ £ 39,656✓ £ 42,722✓ £ 45,787✓
20% £ 75,257✓ £ 79,344✓ £ 83,432✓ £ 87,520✓
25% £ 113,374✓ £ 118,483✓ £ 123,592✓ £ 128,702✓
30% £ 151,112✓ £ 157,248✓ £ 163,384✓ £ 169,520✓
── Participation Incentive Sensitivity (30%, SSEN £250/MWh) ──
£48/HH/yr Kaluza incentive (base £48/HH): P&L £151,112 ✓
£60/HH/yr Kaluza incentive (upper £60/HH): P&L £113,720 ✓
── Thermal Cap Sensitivity (30%, SSEN £250/MWh) ──
Base (10 kWh/day): 1,900 MWh/yr → P&L £151,112 ✓
Sensitivity (15 kWh/day): 2,594 MWh/yr → P&L £341,962 ✓
── Tornado — 30% Penetration | Base profit: £151,112 ──
▲ SSEN £400/MWh (upside) £+285,000
▼ SSEN £150/MWh (downside) £-190,000
▼ Wind farmer fee £0/MWh £-47,500
▼ Incentive £60/HH Kaluza (upper) £-37,392
▲ Wind farmer fee £40/MWh £+28,500
▲ Government subsidy 75% £+18,408
▲ Government subsidy 25% £+6,136
Sandbox¶
In [10]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# ── Load data ─────────────────────────────────────────────────────────────────
df_curt = pd.read_csv('data/curtailment_hourly.csv',
parse_dates=['Timestamp'], index_col='Timestamp')
df_demand = pd.read_csv('data/Residential_demand.csv',
parse_dates=['Timestamp'], index_col='Timestamp')
# ── Hourly means ──────────────────────────────────────────────────────────────
curt_hourly = df_curt['Curtailment_kw'].groupby(df_curt.index.hour).mean() / 1000
demand_hourly = df_demand['Demand_mean_kw'].groupby(df_demand.index.hour).mean()
# ── Normalise to 0-1 ──────────────────────────────────────────────────────────
curt_norm = (curt_hourly - curt_hourly.min()) / (curt_hourly.max() - curt_hourly.min())
demand_norm = (demand_hourly - demand_hourly.min()) / (demand_hourly.max() - demand_hourly.min())
# ── Plot ──────────────────────────────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(12, 5))
hours = np.arange(24)
hour_labels = [f'{h:02d}:00' for h in hours]
ax.plot(hours, curt_norm, color=C2, linewidth=2,
label='Curtailment intensity', zorder=3)
ax.fill_between(hours, curt_norm, alpha=0.15, color=C2)
ax.plot(hours, demand_norm, color=C3, linewidth=2,
label='Residential demand', zorder=3)
ax.fill_between(hours, demand_norm, alpha=0.15, color=C3)
# Callout lines
ax.axvline(13, color=C2, linewidth=1.2, linestyle='--', alpha=0.7)
ax.axvline(17.5, color=C3, linewidth=1.2, linestyle='--', alpha=0.7)
# Annotations
ax.annotate('Curtailment peaks\n13:00 — 78.1 MW mean',
xy=(13, curt_norm[13]), xytext=(10, 0.92),
fontsize=9, color=C2, fontweight='bold',
arrowprops=dict(arrowstyle='->', color=C2, lw=1.2))
ax.annotate('Demand peaks\n17:00–18:00',
xy=(17, demand_norm[17]), xytext=(19, 0.92),
fontsize=9, color=C3, fontweight='bold',
arrowprops=dict(arrowstyle='->', color=C3, lw=1.2))
# Shade the 4-hour window
ax.axvspan(13, 17, alpha=0.08, color=C4, zorder=0)
ax.text(15, 0.05, '4-hour\nwindow', ha='center', fontsize=9,
color=C4, fontweight='bold')
ax.set_xticks(hours[::2])
ax.set_xticklabels(hour_labels[::2], fontsize=9, rotation=45)
ax.set_xlim(0, 23)
ax.set_ylim(0, 1.05)
ax.set_xlabel('Hour of day')
ax.set_ylabel('Normalised intensity')
ax.set_title('Curtailment peaks midday — demand peaks evening\n'
'The 4-hour window Kaluza turns into stored heat',
fontsize=12, fontweight='bold')
ax.legend(fontsize=9, loc='upper left')
plt.tight_layout()
plt.savefig('fig_presentation_slide3_mismatch.png', dpi=150, bbox_inches='tight')
plt.show()
In [14]:
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import numpy as np
penetrations = ['10%', '20%', '30%']
profits = [-2901, 75257, 151112]
colors = [C2, C3, C1]
fig, ax = plt.subplots(figsize=(6, 6))
bars = ax.bar(penetrations, [p/1000 for p in profits],
color=colors, width=0.6, alpha=0.9)
ax.axhline(0, color='#444444', linewidth=0.8, linestyle='--')
for bar, profit in zip(bars, profits):
label = f'£{profit/1000:.0f}k' if profit >= 0 else f'-£{abs(profit)/1000:.0f}k'
ypos = bar.get_height() + 3 if profit >= 0 else bar.get_height() - 8
ax.text(bar.get_x() + bar.get_width()/2, ypos,
label, ha='center', fontsize=13, fontweight='bold',
color='#111111')
hh_labels = ['1,039 HH', '2,077 HH', '3,116 HH']
for i, (bar, label) in enumerate(zip(bars, hh_labels)):
ax.text(bar.get_x() + bar.get_width()/2, -28,
label, ha='center', fontsize=9,
color='#555555')
ax.annotate('Break-even: 15%\n1,558 households',
xy=(0.5, 0), xytext=(0.5, 55),
xycoords=('data', 'data'),
fontsize=9, color='#444444',
ha='center',
arrowprops=dict(arrowstyle='->', color='#888888', lw=0.8))
ax.set_title('Annual accounting profit by Sc.A penetration\n'
'(SSEN £250/MWh, wind farmer fee £25/MWh, 0% subsidy)',
fontsize=11, fontweight='bold', pad=12)
ax.set_xlabel('Sc.A penetration (% of total Orkney households)', fontsize=10)
ax.set_ylabel('Annual accounting profit (£k)', fontsize=10)
ax.set_ylim(-50, 185)
ax.yaxis.set_major_formatter(mticker.FuncFormatter(
lambda x, _: f'£{x:.0f}k'))
plt.tight_layout()
plt.savefig('fig_slide6_profit_bars.png', dpi=150, bbox_inches='tight')
plt.show()
In [15]:
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
RAMP = {2018: 0.30, 2019: 0.70, 2020: 1.00, 2021: 1.00, 2022: 1.00}
ADDITIONS = {2018: 0.30, 2019: 0.40, 2020: 0.30, 2021: 0.00, 2022: 0.00}
YEARS = [2018, 2019, 2020, 2021, 2022]
SSEN_BASE = 250
GEN_FEE = 25
OPEX_PER_HH = 79
FIXED_OPEX = 85500
FUEL_POV = 0.58
CAPEX_A = 225
ONE_OFF = 60
DEVICE_LIFE = 12
scenarios = [
('10%', 1039, '#e8a09a'),
('20%', 2077, '#1D9E75'),
('30%', 3116, '#0F3D2E'),
]
fig, ax = plt.subplots(figsize=(8, 5))
for label, hh, color in scenarios:
n_nonelig = hh - int(hh * FUEL_POV)
total_capex = n_nonelig * CAPEX_A + hh * ONE_OFF
cum = 0.0
cum_vals = []
for yr in YEARS:
ramp = RAMP[yr]
add = ADDITIONS[yr]
rev = hh * ramp * (646/1039) * (SSEN_BASE + GEN_FEE)
if label == '20%':
rev = hh * ramp * (1277/2077) * (SSEN_BASE + GEN_FEE)
elif label == '30%':
rev = hh * ramp * (1900/3116) * (SSEN_BASE + GEN_FEE)
opex = (hh * OPEX_PER_HH + FIXED_OPEX) * ramp
cap = total_capex * add
cum += rev - opex - cap
cum_vals.append(cum / 1000)
ax.plot([str(y) for y in YEARS], cum_vals,
marker='o', label=label, color=color, linewidth=2.5)
ax.axhline(0, color='#444444', linewidth=0.8, linestyle='--')
ax.set_xlabel('Year', fontsize=10)
ax.set_ylabel('Cumulative cash flow (£k)', fontsize=10)
ax.set_title('5-Year Cumulative Cash Flow by Sc.A Penetration\n'
'(SSEN £250/MWh, wind farmer fee £25/MWh, 0% subsidy)',
fontsize=11, fontweight='bold')
ax.legend(title='Penetration', fontsize=9)
ax.yaxis.set_major_formatter(mticker.FuncFormatter(
lambda x, _: f'£{x:,.0f}k'))
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.tight_layout()
plt.savefig('fig_slide6_jcurve.png', dpi=150, bbox_inches='tight')
plt.show()
print("Saved: fig_slide6_jcurve.png")
Saved: fig_slide6_jcurve.png
In [18]:
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import numpy as np
penetrations = ['10%', '20%', '30%']
profits = [-2901, 75257, 151112]
colors = ['#e8a09a', '#1D9E75', '#0F3D2E']
fig, ax = plt.subplots(figsize=(6, 6))
bars = ax.bar(penetrations, [p/1000 for p in profits],
color=colors, width=0.6)
ax.axhline(0, color='#444444', linewidth=0.8, linestyle='--')
for bar, profit in zip(bars, profits):
label = f'£{profit/1000:.0f}k' if profit >= 0 else f'-£{abs(profit)/1000:.0f}k'
ypos = bar.get_height() + 3 if profit >= 0 else bar.get_height() - 8
ax.text(bar.get_x() + bar.get_width()/2, ypos,
label, ha='center', fontsize=13, fontweight='bold',
color='#111111')
hh_labels = ['1,039 HH', '2,077 HH', '3,116 HH']
for bar, label in zip(bars, hh_labels):
ax.text(bar.get_x() + bar.get_width()/2, -28,
label, ha='center', fontsize=9,
color='#555555')
ax.annotate('Break-even: 15%\n1,558 households',
xy=(0.5, 0), xytext=(0.5, 55),
xycoords=('data', 'data'),
fontsize=9, color='#444444',
ha='center',
arrowprops=dict(arrowstyle='->', color='#888888', lw=0.8))
ax.set_title('Annual accounting profit by Sc.A penetration\n'
'(SSEN £250/MWh, wind farmer fee £25/MWh, 0% subsidy)',
fontsize=11, fontweight='bold', pad=12)
ax.set_xlabel('Sc.A penetration (% of total Orkney households)', fontsize=10)
ax.set_ylabel('Annual accounting profit (£k)', fontsize=10)
ax.set_ylim(-50, 185)
ax.yaxis.set_major_formatter(mticker.FuncFormatter(
lambda x, _: f'£{x:.0f}k'))
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.tight_layout()
plt.savefig('fig_slide6_profit_bars.png', dpi=150, bbox_inches='tight')
plt.show()
