Zambia Model

PyPSA Modelling Exercise — Zambian Power System

Overview

In this exercise you will build a stochastic optimisation model using PyPSA.

---

Provided Files

File	Description
Zambia Example.pdf	System diagram with model topology and all static component data
Zambia Data.xlsx	Time-series data (e.g. load profiles, renewable capacity factors)
network_notebook.ipynb	Template notebook — complete your work here

---

Your Task

Using the template notebook, build the PyPSA network.

lesson_number = 5
print(f'lesson{lesson_number}')

lesson5

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=1My8j2qRcjjhVhC5bL657oYTk7-OKDxkE"


else:
    print("Skipping file download.")

import pypsa
import pandas as pd
pypsa.options.api.new_components_api = True
pd.options.plotting.backend = "plotly"

zm = pypsa.Network('zambia_model.xlsx')

# zm.optimize.add_load_shedding()
# zm.optimize.add_load_shedding(p_nom = 10000 , marginal_cost=1000)

zm.optimize()

df_demand = pd.read_excel(
    'Zambia Data.xlsx',
    sheet_name='Demand Data',
    index_col=0,       # first column as index
    header=[0, 1]      # first two rows form a MultiIndex column
)

marginal_cost_dict = {
    'Zambia Coal':{'high':40,'medium':30,'low':20},
    'ZESCO':{'high':45,'medium':35,'low':25}
    }

installed_capacity_dict = {
    'Zambia Coal':{'high':500,'medium':400,'low':300},
    'ZESCO':{'high':2400,'medium':2300,'low':2200}
    }


zm.generators.static.loc[['Zambia Coal','ZESCO'],'marginal_cost']

# ── 1. Define scenarios & probabilities ─────────────────────────────────────
# set_scenarios() modifies the network IN PLACE — do NOT assign its return value
zm.set_scenarios({"Low": 0.2, "Medium": 0.5, "High": 0.3})

# ── 2. Static: marginal costs (Zambia Coal) ──────────────────────────────────
zm.generators.static.loc[("Low",    "Zambia Coal"), "marginal_cost"] = marginal_cost_dict["Zambia Coal"]["low"]
zm.generators.static.loc[("Medium", "Zambia Coal"), "marginal_cost"] = marginal_cost_dict["Zambia Coal"]["medium"]
zm.generators.static.loc[("High",   "Zambia Coal"), "marginal_cost"] = marginal_cost_dict["Zambia Coal"]["high"]

# ── 3. Static: marginal costs (ZESCO) ────────────────────────────────────────
zm.generators.static.loc[("Low",    "ZESCO"), "marginal_cost"] = marginal_cost_dict["ZESCO"]["low"]
zm.generators.static.loc[("Medium", "ZESCO"), "marginal_cost"] = marginal_cost_dict["ZESCO"]["medium"]
zm.generators.static.loc[("High",   "ZESCO"), "marginal_cost"] = marginal_cost_dict["ZESCO"]["high"]

# ── 4. Static: installed capacity / p_nom (Zambia Coal) ──────────────────────
zm.generators.static.loc[("Low",    "Zambia Coal"), "p_nom"] = installed_capacity_dict["Zambia Coal"]["low"]
zm.generators.static.loc[("Medium", "Zambia Coal"), "p_nom"] = installed_capacity_dict["Zambia Coal"]["medium"]
zm.generators.static.loc[("High",   "Zambia Coal"), "p_nom"] = installed_capacity_dict["Zambia Coal"]["high"]

# ── 5. Dynamic: demand load profile (Zambia Load) ────────────────────────────
zm.loads.dynamic.p_set.loc[:, ("Low",    "Zambia Load")] = df_demand[("Low",    "Zambia")].values
zm.loads.dynamic.p_set.loc[:, ("Medium", "Zambia Load")] = df_demand[("Medium", "Zambia")].values
zm.loads.dynamic.p_set.loc[:, ("High",   "Zambia Load")] = df_demand[("High",   "Zambia")].values


# ── 6. Disable cyclic SOC for stochastic solve ──────────────────
# PyPSA bug: cyclic_state_of_charge=True conflicts with scenario MultiIndex indexing.
# Disable it per-scenario so the constraint builder doesn't hit a 2D boolean issue.
zm.storage_units.static.loc[("Low",    "New Zambia HPP"), "cyclic_state_of_charge"] = False
zm.storage_units.static.loc[("Medium", "New Zambia HPP"), "cyclic_state_of_charge"] = False
zm.storage_units.static.loc[("High",   "New Zambia HPP"), "cyclic_state_of_charge"] = False

# ── 14. Solve ─────────────────────────────────────────────────────────────────
zm.optimize(solver_name="highs")

zm.buses.static

zm.generators.static.p_nom_opt

zm.generators.dynamic.p_max_pu['Low'].plot()

zm.loads.dynamic.p_set['Low'].plot()
zm.loads.dynamic.p_set['Medium'].plot()
zm.loads.dynamic.p_set['High'].plot()

import plotly.graph_objects as go

fig = go.Figure()
for scenario in ['Low', 'Medium', 'High']:
    s = zm.loads.dynamic.p_set[scenario].sum(axis=1)
    fig.add_trace(go.Scatter(x=s.index, y=s.values, name=scenario))

fig.update_layout(xaxis_title='Time', yaxis_title='Load (MW)')
fig.show()

zm.generators.dynamic.p.sum()

# ── DIAGNOSTIC: Why is p_nom_opt the same across scenarios? ──────────────────

print("=" * 60)
print("1. EXTENDABLE vs FIXED generators")
print("=" * 60)
# p_nom_extendable=True → investment decision (first-stage) → MUST be same across scenarios
# p_nom_extendable=False → p_nom_opt = p_nom (just reflects the input)
ext = zm.generators.static[["p_nom", "p_nom_extendable", "p_nom_opt"]].xs("Low", level="scenario")
print(ext.to_string())

print()
print("=" * 60)
print("2. p_nom INPUT per scenario (was it set differently?)")
print("=" * 60)
# Shows which generators actually received different p_nom per scenario
p_nom_by_scenario = zm.generators.static["p_nom"].unstack("scenario")[["Low", "Medium", "High"]]
print(p_nom_by_scenario.to_string())

print()
print("=" * 60)
print("3. p_nom_opt OUTPUT per scenario")
print("=" * 60)
p_nom_opt_by_scenario = zm.generators.static["p_nom_opt"].unstack("scenario")[["Low", "Medium", "High"]]
print(p_nom_opt_by_scenario.to_string())

print()
print("=" * 60)
print("4. DISPATCH (p) — this SHOULD differ across scenarios")
print("   Mean dispatch per generator per scenario:")
print("=" * 60)
for scenario in ["Low", "Medium", "High"]:
    print(f"\n  Scenario: {scenario}")
    print(zm.generators.dynamic.p[scenario].mean().to_string())

import plotly.graph_objects as go

gen_total = zm.generators.dynamic.p.sum().unstack("scenario")[["Low", "Medium", "High"]]

fig = go.Figure()
for scenario in ["Low", "Medium", "High"]:
    fig.add_trace(go.Bar(
        name=scenario,
        x=gen_total.index,
        y=gen_total[scenario]
    ))

fig.update_layout(
    barmode="group",
    title="Total Generation by Scenario",
    xaxis_title="Generator",
    yaxis_title="Total Generation (MWh)",
    xaxis_tickangle=-45
)
fig.show()