Exercise 2: Capacity Expansion Planning#
Prepare Google Colab Environment#
#@title Connect to Google Drive {display-mode:"form"}
CONNECT_TO_DRIVE = False #@param {type:"boolean"}
import os
if CONNECT_TO_DRIVE:
from google.colab import drive
# Mount Google Drive
drive.mount('/content/drive')
# Define the desired working directory path
working_dir = '/content/drive/MyDrive/hello-pypsa'
# Create the directory if it doesn't exist
if not os.path.exists(working_dir):
os.makedirs(working_dir)
print(f"Directory '{working_dir}' created.")
else:
print(f"Directory '{working_dir}' already exists.")
# Change the current working directory
os.chdir(working_dir)
print(f"Current working directory: {os.getcwd()}")
else:
print("Not connecting to Google Drive.")
import os
#@title Install Packages {display-mode:"form"}
INSTALL_PACKAGES = False #@param {type:"boolean"}
# Check if packages have already been installed in this session to prevent re-installation
if INSTALL_PACKAGES and not os.environ.get('PYPSA_PACKAGES_INSTALLED'):
!pip install pypsa pypsa[excel] folium mapclassify cartopy
!pip install git+https://github.com/PriyeshGosai/pypsa_network_viewer.git
os.environ['PYPSA_PACKAGES_INSTALLED'] = 'true'
elif not INSTALL_PACKAGES:
print("Skipping package installation.")
else:
print("PyPSA packages are already installed for this session.")
#@title Download the file for this notebook {display-mode:"form"}
DOWNLOAD_FILE = False #@param {type:"boolean"}
if DOWNLOAD_FILE:
!gdown "https://drive.google.com/uc?id=14DYNjwrDy3XFIzba6Cv2C9hsk6LdZsuw"
!gdown "https://drive.google.com/uc?id=1ke-wMEQZCpf8erPGErT0lLdGce7EocJn"
else:
print("Skipping file download.")
Exercise 2#
print("PyPSA Model")
PyPSA Model
import pypsa
import geopandas as gpd
pypsa.options.api.new_components_api = True
import matplotlib.pyplot as plt
import pandas as pd
pd.options.plotting.backend = 'plotly'
gpkg_file_name = "zaf_admin1.gpkg"
network_file = 'za_model.xlsx'
gdf_index = 'adm1_name'
n_days = 10
gdf = gpd.read_file(gpkg_file_name)
# gdf.loc['adm1_name', '7'] = 'Northern Cape'
n = pypsa.Network(network_file)
gdf['name'] = gdf['adm1_name'].str.lower().str.replace(r'[ -]', '_', regex=True)
gdf.set_index('name',inplace=True)
n.add("Shape",
gdf.index,
geometry=gdf['geometry'],
component='Bus',
idx=gdf['adm1_name'])
n.generators.static.capital_cost *= pypsa.costs.annuity(0.08, 20)
if n_days is not None:
n_hours = n_days * 24
n.generators.static.capital_cost *= pypsa.costs.annuity(0.08, 20)*n_days/365
n.set_snapshots(n.snapshots[:n_hours])
n.global_constraints.static.loc['carbon_constraint','constant'] *= n_days/365
n.optimize(include_objective_constant=True)
import plotly.express as px
df = n.generators.static[['p_nom', 'p_nom_opt','p_nom_min']].reset_index()
fig = px.bar(
df.melt(id_vars='name', value_vars=['p_nom', 'p_nom_opt']),
x='name',
y='value',
color='variable',
barmode='group',
title='Generator Capacity',
labels={'value': 'Power (MW)', 'variable': 'Capacity Type', 'name': 'Generator'}
)
fig.update_layout(xaxis_tickangle=-45)
fig.show()
n.export_to_netcdf('reults.nc')
n.generators.dynamic.p.plot.area(stacked = True)
from matplotlib.patches import Wedge
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
# --- Palette ---
BG_COLOR = "#ffffff" # dark navy background
LAND_COLOR = "#cecece" # slightly lighter navy for land
BORDER_COLOR = '#0f3460' # subtle blue borders
LINK_COLOR = '#e94560' # coral/red for HVDC links
LINE_COLOR = '#4cc9f0' # light blue for AC lines
CARRIER_COLORS = {
'coal': '#e76f51', # burnt orange
'gas': '#f4a261', # warm amber
'diesel': '#e9c46a', # yellow
'solar': '#ffd166', # bright yellow
'wind': '#06d6a0', # teal green
'nuclear': '#8338ec', # purple
'hydro': '#4cc9f0', # sky blue
'load shedding': '#ef233c', # red
}
def get_carrier_color(carrier):
"""Fallback to a muted grey for unknown carriers."""
return CARRIER_COLORS.get(carrier, '#adb5bd')
# --- Figure ---
fig, ax = plt.subplots(figsize=(12, 12), facecolor=BG_COLOR)
ax.set_facecolor(BG_COLOR)
# 1. Background shapes
n.shapes.static.plot(ax=ax, color=LAND_COLOR, edgecolor=BORDER_COLOR, linewidth=0.7)
# 2. Links (HVDC)
for name, row in n.links.static.iterrows():
bus0 = n.buses.static.loc[row['bus0']]
bus1 = n.buses.static.loc[row['bus1']]
ax.plot(
[bus0['x'], bus1['x']],
[bus0['y'], bus1['y']],
linewidth=np.log1p(row['p_nom']) * 0.5,
color=LINK_COLOR,
alpha=0.8,
zorder=2,
linestyle='--'
)
# 3. Lines (AC)
for name, row in n.lines.static.iterrows():
bus0 = n.buses.static.loc[row['bus0']]
bus1 = n.buses.static.loc[row['bus1']]
ax.plot(
[bus0['x'], bus1['x']],
[bus0['y'], bus1['y']],
linewidth=np.log1p(row['p_nom']) * 0.5,
color=LINE_COLOR,
alpha=0.5,
zorder=3,
)
# 4. Pie charts
grouped = n.generators.static.groupby(['bus', 'carrier'])['p_nom'].sum()
scale = 0.5
for bus, row in n.buses.static.iterrows():
if bus not in grouped.index.get_level_values('bus'):
continue
data = grouped[bus]
total = data.sum()
radius = np.log1p(total / 1e3) * scale
start = 0
for carrier, val in data.items():
angle = 360 * val / total
wedge = Wedge(
(row['x'], row['y']), radius,
start, start + angle,
facecolor=get_carrier_color(carrier), # changed from color=
alpha=0.95,
zorder=5,
linewidth=0.3,
edgecolor=BG_COLOR
)
ax.add_patch(wedge)
start += angle
ax.annotate(
f'{int(total)} MW',
(row['x'], row['y'] + radius + 0.05),
fontsize=7,
ha='center',
color='black',
zorder=6
)
# 5. Legend
carriers_present = n.generators.static['carrier'].unique()
legend_patches = [
mpatches.Patch(color=get_carrier_color(c), label=c)
for c in carriers_present
]
legend_patches += [
plt.Line2D([0], [0], color=LINE_COLOR, linewidth=2, linestyle='-', label='AC line'),
plt.Line2D([0], [0], color=LINK_COLOR, linewidth=2, linestyle='--', label='HVDC link'),
]
ax.legend(
handles=legend_patches,
loc='lower right',
framealpha=0.2,
facecolor=BG_COLOR,
edgecolor=BORDER_COLOR,
labelcolor='black',
fontsize=9
)
ax.set_title('Generator capacity by bus and carrier', color='black', fontsize=14, pad=12)
ax.set_axis_off()
plt.tight_layout()
plt.show()
from matplotlib.patches import Wedge
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
# --- Palette ---
BG_COLOR = "#ffffff" # dark navy background
LAND_COLOR = "#cecece" # slightly lighter navy for land
BORDER_COLOR = '#0f3460' # subtle blue borders
LINK_COLOR = '#e94560' # coral/red for HVDC links
LINE_COLOR = '#4cc9f0' # light blue for AC lines
CARRIER_COLORS = {
'coal': '#e76f51', # burnt orange
'gas': '#f4a261', # warm amber
'diesel': '#e9c46a', # yellow
'solar': '#ffd166', # bright yellow
'wind': '#06d6a0', # teal green
'nuclear': '#8338ec', # purple
'hydro': '#4cc9f0', # sky blue
'load shedding': '#ef233c', # red
}
def get_carrier_color(carrier):
"""Fallback to a muted grey for unknown carriers."""
return CARRIER_COLORS.get(carrier, '#adb5bd')
# --- Figure ---
fig, ax = plt.subplots(figsize=(12, 12), facecolor=BG_COLOR)
ax.set_facecolor(BG_COLOR)
# 1. Background shapes
n.shapes.static.plot(ax=ax, color=LAND_COLOR, edgecolor=BORDER_COLOR, linewidth=0.7)
# 2. Links (HVDC)
for name, row in n.links.static.iterrows():
bus0 = n.buses.static.loc[row['bus0']]
bus1 = n.buses.static.loc[row['bus1']]
ax.plot(
[bus0['x'], bus1['x']],
[bus0['y'], bus1['y']],
linewidth=np.log1p(row['p_nom']) * 0.5,
color=LINK_COLOR,
alpha=0.8,
zorder=2,
linestyle='--'
)
# 3. Lines (AC)
for name, row in n.lines.static.iterrows():
bus0 = n.buses.static.loc[row['bus0']]
bus1 = n.buses.static.loc[row['bus1']]
ax.plot(
[bus0['x'], bus1['x']],
[bus0['y'], bus1['y']],
linewidth=np.log1p(row['p_nom']) * 0.5,
color=LINE_COLOR,
alpha=0.5,
zorder=3,
)
# 4. Pie charts
grouped = n.generators.static.groupby(['bus', 'carrier'])['p_nom'].sum()
scale = 0.5
for bus, row in n.buses.static.iterrows():
if bus not in grouped.index.get_level_values('bus'):
continue
data = grouped[bus]
total = data.sum()
radius = np.log1p(total / 1e3) * scale
start = 0
for carrier, val in data.items():
angle = 360 * val / total
wedge = Wedge(
(row['x'], row['y']), radius,
start, start + angle,
facecolor=get_carrier_color(carrier), # changed from color=
alpha=0.95,
zorder=5,
linewidth=0.3,
edgecolor=BG_COLOR
)
ax.add_patch(wedge)
start += angle
ax.annotate(
f'{int(total)} MW',
(row['x'], row['y'] + radius + 0.05),
fontsize=7,
ha='center',
color='black',
zorder=6
)
# 5. Legend
carriers_present = n.generators.static['carrier'].unique()
legend_patches = [
mpatches.Patch(color=get_carrier_color(c), label=c)
for c in carriers_present
]
legend_patches += [
plt.Line2D([0], [0], color=LINE_COLOR, linewidth=2, linestyle='-', label='AC line'),
plt.Line2D([0], [0], color=LINK_COLOR, linewidth=2, linestyle='--', label='HVDC link'),
]
ax.legend(
handles=legend_patches,
loc='lower right',
framealpha=0.2,
facecolor=BG_COLOR,
edgecolor=BORDER_COLOR,
labelcolor='black',
fontsize=9
)
ax.set_title('Generator capacity by bus and carrier', color='black', fontsize=14, pad=12)
ax.set_axis_off()
plt.tight_layout()
plt.show()
from matplotlib.patches import Circle
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
# --- Palette ---
BG_COLOR = "#ffffff"
LAND_COLOR = "#cecece"
BORDER_COLOR = '#0f3460'
LOAD_COLOR = "#4558e9"
LINE_COLOR = '#4cc9f0'
# --- Data: map bus name -> annual load ---
load_totals = n.loads.dynamic.p_set.sum()
bus_map = n.loads.static['bus'] # load name -> bus name
load_by_bus = {bus_map[load]: val for load, val in load_totals.items()}
print("load_by_bus:", load_by_bus) # sanity check — remove once confirmed working
# --- Figure ---
fig, ax = plt.subplots(figsize=(12, 12), facecolor=BG_COLOR)
ax.set_facecolor(BG_COLOR)
# 1. Background shapes
n.shapes.static.plot(ax=ax, color=LAND_COLOR, edgecolor=BORDER_COLOR, linewidth=0.7)
# 2. Lines (AC)
for name, row in n.lines.static.iterrows():
bus0 = n.buses.static.loc[row['bus0']]
bus1 = n.buses.static.loc[row['bus1']]
ax.plot(
[bus0['x'], bus1['x']],
[bus0['y'], bus1['y']],
linewidth=np.log1p(row['p_nom']) * 0.5,
color=LINE_COLOR,
alpha=0.5,
zorder=3,
)
# 3. Load circles
scale = 0.4
for bus, row in n.buses.static.iterrows():
if bus not in load_by_bus:
continue
total = load_by_bus[bus]
radius = np.log1p(total / 1e6) * scale
circle = Circle(
(row['x'], row['y']), radius,
facecolor=LOAD_COLOR,
alpha=0.75,
zorder=5,
linewidth=0.5,
edgecolor='white'
)
ax.add_patch(circle)
ax.annotate(
f'{total/1e6:.1f} TWh',
(row['x'], row['y'] + radius + 0.05),
fontsize=7,
ha='center',
color='black',
zorder=6
)
# 4. Legend
legend_handles = [
mpatches.Patch(color=LOAD_COLOR, alpha=0.75, label='Annual load'),
plt.Line2D([0], [0], color=LINE_COLOR, linewidth=2, linestyle='-', label='AC line'),
]
ax.legend(
handles=legend_handles,
loc='lower right',
framealpha=0.2,
facecolor=BG_COLOR,
edgecolor=BORDER_COLOR,
labelcolor='black',
fontsize=9
)
ax.set_title('Annual load by province', color='black', fontsize=14, pad=12)
ax.set_axis_off()
plt.tight_layout()
plt.show()
from matplotlib.patches import Wedge
import numpy as np
fig, ax = plt.subplots(figsize=(10, 10))
# 1. Background shapes
n.shapes.static.plot(ax=ax, color='white', edgecolor='grey', linewidth=0.5)
# 2. Links (HVDC etc.)
for name, row in n.links.static.iterrows():
bus0 = n.buses.static.loc[row['bus0']]
bus1 = n.buses.static.loc[row['bus1']]
ax.plot(
[bus0['x'], bus1['x']],
[bus0['y'], bus1['y']],
linewidth=np.log1p(row['p_nom']) * 0.5,
color='darkorange',
alpha=0.6,
zorder=2,
linestyle='--'
)
# 3. Lines (AC)
for name, row in n.lines.static.iterrows():
bus0 = n.buses.static.loc[row['bus0']]
bus1 = n.buses.static.loc[row['bus1']]
ax.plot(
[bus0['x'], bus1['x']],
[bus0['y'], bus1['y']],
linewidth=np.log1p(row['p_nom']) * 0.5,
color='steelblue',
alpha=0.6,
zorder=3,
)
# 4. Pie charts on top
carrier_colors = {c: plt.cm.tab10(i) for i, c in enumerate(n.generators.static['carrier'].unique())}
grouped = n.generators.static.groupby(['bus', 'carrier'])['p_nom'].sum()
scale = 0.5
for bus, row in n.buses.static.iterrows():
if bus not in grouped.index.get_level_values('bus'):
continue
data = grouped[bus]
total = data.sum()
radius = np.log1p(total / 1e3) * scale
start = 0
for carrier, val in data.items():
angle = 360 * val / total
wedge = Wedge(
(row['x'], row['y']), radius,
start, start + angle,
color=carrier_colors[carrier], alpha=1,
zorder=5
)
ax.add_patch(wedge)
start += angle
ax.annotate(
f'{int(total)} MW',
(row['x'], row['y'] + radius + 0.05),
fontsize=7,
ha='center',
zorder=6
)
# Legend
for carrier, color in carrier_colors.items():
ax.plot([], [], 'o', color=color, label=carrier)
ax.legend(loc='lower right')
ax.set_title('Generator capacity by bus and carrier')
ax.set_axis_off()
plt.tight_layout()
plt.show()
eb = (
n.statistics.energy_balance(
groupby=["bus", "carrier"],
components=["Generator", "Load", "StorageUnit"],
)
.groupby(["bus", "carrier"])
.sum()
)
line_flow = n.lines.dynamic.p0.sum(axis=0)
link_flow = n.links.dynamic.p0.sum(axis=0)
# view_state = {}
# view_state["zoom"] = 6
# view_state["pitch"] = 35 # Up/down angle relative to the map's plane
# map = n.explore(
# view_state=view_state,
# map_style="dark",
# bus_size=eb, # MWh -> km²
# bus_split_circle=True,
# bus_size_max=7000, # km²
# line_color="yellow",
# line_width=line_flow, # MWh -> km
# link_width=link_flow, # MWh -> km
# line_flow=line_flow, # MWh -> km
# link_flow=link_flow, # MWh -> km
# branch_width_max=16, # km
# auto_scale=True,
# bus_columns=["v_nom"],
# line_columns=["s_nom"],
# link_columns=["p_nom"],
# arrow_size_factor=2,
# tooltip=True, # disabled here for technical limits of mkdocs-jupyter plugin
# )
n.generators.static.columns
n.generators.static[['p_nom_opt','p_nom']].round(0).astype(int)
n.statistics.system_cost()
n.generators.dynamic.p.plot()
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
utilisation = np.abs(n.links_t.p0 / n.links.static.p_nom)
bins = [0, 0.25, 0.50, 0.75, 1.0]
labels = ['0–25%', '25–50%', '50–75%', '75–100%']
colors = ['#4cc9f0', '#06d6a0', '#f4a261', '#e76f51']
# Build a DataFrame of percentages: rows=links, columns=bands
pct_df = pd.DataFrame(index=utilisation.columns, columns=labels, dtype=float)
for link in utilisation.columns:
counts = pd.cut(utilisation[link], bins=bins, labels=labels, include_lowest=True).value_counts()
pct_df.loc[link] = (counts.reindex(labels) / len(utilisation) * 100).values
# Stacked bar chart
fig, ax = plt.subplots(figsize=(10, 5))
bottom = np.zeros(len(pct_df))
for label, color in zip(labels, colors):
bars = ax.bar(pct_df.index, pct_df[label], bottom=bottom, label=label,
color=color, edgecolor='white', linewidth=0.5)
for bar, val in zip(bars, pct_df[label]):
if val > 3:
ax.text(bar.get_x() + bar.get_width() / 2,
bar.get_y() + bar.get_height() / 2,
f'{val:.1f}%', ha='center', va='center', fontsize=8, color='black')
bottom += pct_df[label].values
ax.set_ylabel('% of snapshots')
ax.set_xlabel('Link')
ax.set_ylim(0, 100)
ax.set_title('Link utilisation distribution')
ax.legend(title='Utilisation band', bbox_to_anchor=(1.01, 1), loc='upper left')
ax.tick_params(axis='x', rotation=30)
plt.tight_layout()
plt.show()
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
utilisation = n.links.dynamic.p0 / n.links.static.p_nom
bins = [-1.0, -0.75, -0.50, -0.25, 0.0, 0.25, 0.50, 0.75, 1.0]
labels = ['-100–-75%', '-75–-50%', '-50–-25%', '-25–0%',
'0–25%', '25–50%', '50–75%', '75–100%']
# Reds for reverse (dark → light approaching zero), blues for forward (light → dark)
colors = ['#7b0000', '#c1121f', '#e76f51', '#ffb4a2', # reverse: dark to light
'#a8dadc', '#4cc9f0', '#1d7ec1', '#03045e'] # forward: light to dark
pct_df = pd.DataFrame(index=utilisation.columns, columns=labels, dtype=float)
for link in utilisation.columns:
counts = pd.cut(utilisation[link], bins=bins, labels=labels, include_lowest=True).value_counts()
pct_df.loc[link] = (counts.reindex(labels) / len(utilisation) * 100).values
fig, ax = plt.subplots(figsize=(12, 8))
bottom = np.zeros(len(pct_df))
for label, color in zip(labels, colors):
bars = ax.bar(pct_df.index, pct_df[label], bottom=bottom, label=label,
color=color, edgecolor='white', linewidth=0.5)
for bar, val in zip(bars, pct_df[label]):
if val > 3:
ax.text(bar.get_x() + bar.get_width() / 2,
bar.get_y() + bar.get_height() / 2,
f'{val:.1f}%', ha='center', va='center', fontsize=8, color='white')
bottom += pct_df[label].values
ax.set_ylabel('% of snapshots')
ax.set_xlabel('Link')
ax.set_ylim(0, 100)
ax.set_title('Link utilisation distribution')
ax.legend(title='Utilisation band', bbox_to_anchor=(1.01, 1), loc='upper left')
ax.tick_params(axis='x', rotation=90)
plt.tight_layout()
plt.show()
import plotly.graph_objects as go
import numpy as np
import pandas as pd
utilisation = n.links.dynamic.p0 / n.links.static.p_nom
bins = [-1.0, -0.75, -0.50, -0.25, 0.0, 0.25, 0.50, 0.75, 1.0]
labels = ['-100–-75%', '-75–-50%', '-50–-25%', '-25–0%',
'0–25%', '25–50%', '50–75%', '75–100%']
colors = ['#7b0000', '#c1121f', '#e76f51', '#ffb4a2',
'#a8dadc', '#4cc9f0', '#1d7ec1', '#03045e']
pct_df = pd.DataFrame(index=utilisation.columns, columns=labels, dtype=float)
for link in utilisation.columns:
counts = pd.cut(utilisation[link], bins=bins, labels=labels, include_lowest=True).value_counts()
pct_df.loc[link] = (counts.reindex(labels) / len(utilisation) * 100).values
fig = go.Figure()
for label, color in zip(labels, colors):
fig.add_trace(go.Bar(
name=label,
x=pct_df.index,
y=pct_df[label],
marker_color=color,
marker_line_color='white',
marker_line_width=0.5,
text=[f'{v:.1f}%' if v > 3 else '' for v in pct_df[label]],
textposition='inside',
textfont=dict(color='white', size=11),
hovertemplate='<b>%{x}</b><br>' + label + ': %{y:.1f}%<extra></extra>',
))
fig.update_layout(
barmode='stack',
title='Link utilisation distribution',
xaxis_title='Link',
yaxis_title='% of snapshots',
yaxis=dict(range=[0, 100]),
legend_title='Utilisation band',
legend=dict(x=1.01, y=1, xanchor='left'),
xaxis_tickangle=-90,
plot_bgcolor='white',
height=550,
)
fig.show()
(n.links.dynamic.p0/n.links.static.p_nom).describe()
# from pypsa_network_viewer import html_network , generate_template
# html_file = html_network(
# n,
# file_name='exercise_2_network.html',
# title='Exercise 2 Analysis'
# )