Tutorial 3: Capacity expansion modelling#

Introduction#

This tutorial demonstrates how to perform capacity expansion modeling using PyPSA (Python for Power System Analysis). We’ll build a simple power system model that optimizes the deployment of renewable energy sources (wind and solar) along with storage capacity to meet electricity demand at minimum cost.

Data preparation#

First, we need to prepare our input data. We’ll process load data and renewable generation profiles.

Step 1: Load data processing#

# imports and settings for the plotting
from pypsa import Network
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

sns.set_theme("notebook", style="white")
plt.rc("axes.spines", top=False, right=False)
# Load the data
RAW_FILE = "data/merit_order.csv"
OUTPUT_FILE = "data/load.csv"
NORMALIZE = False
RESET_INDEX_YEAR = 2019

data = pd.read_csv(RAW_FILE, index_col=0, header=0)
data = data[[col for col in data.columns if ".input" in col]].sum(axis=1)

if NORMALIZE:
    data = data / data.max()

if RESET_INDEX_YEAR:
    data.index = pd.date_range(
        start=f"{RESET_INDEX_YEAR}-01-01", periods=len(data), freq="h"
    )

data = data.to_frame()
data.index.name = "time"
data.columns = ["load"]
data.to_csv(OUTPUT_FILE)

Step 2: Combining timeseries data#

We combine the timeseries to an uniform data format for ease of use in our optimisation model.

# Load the data files

LOAD_FP = "data/load.csv"
SOLAR_FP = "data/solar.csv"
WIND_FP = "data/wind.csv"


load_data = pd.read_csv(
    LOAD_FP, index_col=0, parse_dates=True, header=0, names=["load"]
)
solar_data = pd.read_csv(SOLAR_FP, header=0)
wind_data = pd.read_csv(WIND_FP, header=0)

# Combine based on load timestamps
combined_data = pd.DataFrame(index=load_data.index)
combined_data.index.name = "time"
combined_data["load"] = load_data["load"]
combined_data["solar"] = solar_data["electricity"].values
combined_data["wind"] = wind_data["electricity"].values

# Save the combined timeseries
combined_data.to_csv("data/timeseries.csv")

ts = combined_data

Network creation and optimization#

Now we’ll create our network model. We’ll analyze three different scenarios with varying cost ratios between generation and storage. We will first define some functions that reuse.

def create_network(
    ts,
    costs,
    generation_cost_multiplier=1.0,
    storage_cost_multiplier=1.0,
    add_flow_battery=False,
    limit_renewable_energy=False,
) -> Network:

    # Adjust costs based on multipliers
    costs.loc[["wind", "solar"], "capital_cost"] *= generation_cost_multiplier
    costs.loc[
        ["battery inverter", "battery storage"], "capital_cost"
    ] *= storage_cost_multiplier

    # initialize network
    n = Network()
    n.set_snapshots(ts.index)
    n.add("Bus", "electricity")

    # add load
    n.add(
        "Load",
        "demand",
        bus="electricity",
        p_set=ts.load,
    )

    # add solar and wind
    tech_limits = {
        "wind": 150e3,
        "solar": 300e3,
    }
    for tech in ["wind", "solar"]:
        n.add(
            "Generator",
            tech,
            bus="electricity",
            p_max_pu=ts[tech],
            capital_cost=costs.at[tech, "capital_cost"],
            marginal_cost=costs.at[tech, "marginal_cost"],
            p_nom_extendable=True,
            p_nom_max=tech_limits[tech] if limit_renewable_energy else None,
        )

    # add short-duration battery storage
    EP_RATIO = 6
    n.add(
        "StorageUnit",
        "battery storage",
        bus="electricity",
        max_hours=EP_RATIO,
        capital_cost=costs.at["battery inverter", "capital_cost"]
        + EP_RATIO * costs.at["battery storage", "capital_cost"],
        efficiency_store=costs.at["battery inverter", "efficiency"],
        efficiency_dispatch=costs.at["battery inverter", "efficiency"],
        p_nom_extendable=True,
        cyclic_state_of_charge=True,
    )

    if add_flow_battery:
        # add long-duration flow battery
        FLOW_EP_RATIO = 100
        n.add(
            "StorageUnit",
            "flow battery",
            bus="electricity",
            max_hours=FLOW_EP_RATIO,
            capital_cost=costs.at["flow battery power", "capital_cost"]
            + FLOW_EP_RATIO * costs.at["flow battery energy", "capital_cost"],
            efficiency_store=costs.at["flow battery power", "efficiency"],
            efficiency_dispatch=costs.at["flow battery power", "efficiency"],
            p_nom_extendable=True,
            cyclic_state_of_charge=True,
        )

    return n

Baseline scenario analysis#

For the baseline scenario, we assume only lithium ion batteries to be available.

# load the cost data
COSTS_FP = "data/costs.csv"
SOLVER_OPTIONS = {
    "primal_feasibility_tolerance": 1e-5,
    "dual_feasibility_tolerance": 1e-5,
}
# we decrease the tolerance to speed up the optimization (default is 1e-7)

costs = pd.read_csv(COSTS_FP, index_col=0, header=0)

n = create_network(ts, costs, add_flow_battery=False)

n.optimize(
    solver_name="highs",
    solver_options=SOLVER_OPTIONS,
)

Hide code cell output

/home/thesethtruth/.cache/pypoetry/virtualenvs/esi-course-6syH-QiP-py3.12/lib/python3.12/site-packages/pypsa/io.py:1096: FutureWarning: Downcasting object dtype arrays on .fillna, .ffill, .bfill is deprecated and will change in a future version. Call result.infer_objects(copy=False) instead. To opt-in to the future behavior, set `pd.set_option('future.no_silent_downcasting', True)`
  df = df.fillna(attrs["default"].to_dict())
/home/thesethtruth/.cache/pypoetry/virtualenvs/esi-course-6syH-QiP-py3.12/lib/python3.12/site-packages/pypsa/io.py:1096: FutureWarning: Downcasting object dtype arrays on .fillna, .ffill, .bfill is deprecated and will change in a future version. Call result.infer_objects(copy=False) instead. To opt-in to the future behavior, set `pd.set_option('future.no_silent_downcasting', True)`
  df = df.fillna(attrs["default"].to_dict())
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
/home/thesethtruth/.cache/pypoetry/virtualenvs/esi-course-6syH-QiP-py3.12/lib/python3.12/site-packages/linopy/common.py:154: UserWarning: coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.
  warn(
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options:
 - primal_feasibility_tolerance: 1e-05
 - dual_feasibility_tolerance: 1e-05
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 16/16 [00:00<00:00, 41.47it/s]
Writing continuous variables.: 100%|██████████| 6/6 [00:00<00:00, 99.18it/s]
INFO:linopy.io: Writing time: 0.47s
Running HiGHS 1.8.1 (git hash: 4a7f24a): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-03, 6e+00]
  Cost   [5e-01, 2e+03]
  Bound  [0e+00, 0e+00]
  RHS    [2e+04, 1e+05]
Presolving model
56913 rows, 39396 cols, 144459 nonzeros  0s
56913 rows, 39396 cols, 144459 nonzeros  0s
Presolve : Reductions: rows 56913(-48210); columns 39396(-4407); elements 144459(-52617)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 8760(2.88434e+08) 0s
      17392     8.3806783045e+08 Pr: 12032(1.38802e+11); Du: 0(0.000685959) 5s
      36450     1.2025215527e+09 Pr: 0(0); Du: 0(1.21225e-12) 10s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-jr192ado
Model status        : Optimal
Simplex   iterations: 36450
Objective value     :  1.2025215527e+09
Relative P-D gap    :  2.7757174913e-15
HiGHS run time      :         10.36
Writing the solution to /tmp/linopy-solve-5rhrhsfp.sol
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 43803 primals, 105123 duals
Objective: 1.20e+09
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
('ok', 'optimal')

Plot the energy system configuration#

fig, ax = plt.subplots(1, 1, figsize=(10, 3))
plt.subplots_adjust(hspace=0.4)

capacity_data = pd.concat(
    [
        n.generators.p_nom_opt,
        n.storage_units.p_nom_opt,
    ],
    axis=1,
).div(1e3)

capacity_data.columns = ["generation", "storage"]

capacity_data.plot(kind="barh", ax=ax)
ax.set_title("Installed Capacities: Baseline scenario")
ax.set_ylabel("Capacity (GW)")
Text(0, 0.5, 'Capacity (GW)')
_images/1574a3a3ed7fed7c6feafeb1e53379af597112dfc9aa75106ddba4f1ba5ba789.png

The baseline scenario shows very large capacities;

  • 500 GW solar

  • 230 GW wind

  • and 130 GW storage

These values exceed what could physically fit in the Netherlands. This oversizing occurs because our model treats the Netherlands as an isolated system without connections to neighboring countries, requiring full self-sufficiency in power generation.

Baseline with renewable generation limits#

n = create_network(ts, costs, add_flow_battery=False, limit_renewable_energy=True)

n.optimize(
    solver_name="highs",
    solver_options=SOLVER_OPTIONS,
)

Hide code cell output

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
/home/thesethtruth/.cache/pypoetry/virtualenvs/esi-course-6syH-QiP-py3.12/lib/python3.12/site-packages/linopy/common.py:154: UserWarning: coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.
  warn(
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options:
 - primal_feasibility_tolerance: 1e-05
 - dual_feasibility_tolerance: 1e-05
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 16/16 [00:00<00:00, 42.99it/s]
Writing continuous variables.: 100%|██████████| 6/6 [00:00<00:00, 82.10it/s]
INFO:linopy.io: Writing time: 0.48s
Running HiGHS 1.8.1 (git hash: 4a7f24a): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-03, 6e+00]
  Cost   [5e-01, 2e+03]
  Bound  [0e+00, 0e+00]
  RHS    [2e+04, 3e+05]
Presolving model
56913 rows, 39396 cols, 144459 nonzeros  0s
56913 rows, 39396 cols, 144459 nonzeros  0s
Presolve : Reductions: rows 56913(-48212); columns 39396(-4407); elements 144459(-52619)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 8760(2.88434e+08) 0s
      22648     1.0572706864e+09 Pr: 9883(7.01578e+10); Du: 0(0.00171288) 5s
      31363     1.2643315278e+09 Pr: 2006(3.50648e+07); Du: 0(0.00374069) 11s
      39571     1.7644238200e+09 Pr: 0(0); Du: 0(6.7335e-13) 11s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-b6c1or7u
Model status        : Optimal
Simplex   iterations: 39571
Objective value     :  1.7644238200e+09
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 43803 primals, 105125 duals
Objective: 1.76e+09
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
Relative P-D gap    :  2.9051973753e-14
HiGHS run time      :         10.96
Writing the solution to /tmp/linopy-solve-anmbhe_6.sol
('ok', 'optimal')

Plot the energy system configuration#

fig, ax = plt.subplots(1, 1, figsize=(10, 3))
plt.subplots_adjust(hspace=0.4)

capacity_data = pd.concat(
    [
        n.generators.p_nom_opt,
        n.storage_units.p_nom_opt,
    ],
    axis=1,
).div(1e3)

capacity_data.columns = ["generation", "storage"]

capacity_data.plot(kind="barh", ax=ax)
ax.set_title("Installed Capacities: Baseline with limited renewable energy")
ax.set_ylabel("Capacity (GW)")
Text(0, 0.5, 'Capacity (GW)')
_images/dc5148284125c7027c46172d81070d5135b3198f13ceb95b8fe058824b535a81.png

Impact of generation capacity limits#

With upper limits applied to renewable generation, we see:

  • 300 GW solar (at limit)

  • 150 GW (at limit)

  • 480 GW lithium-ion battery storage

The model responds to these generation constraints by heavily investing in storage capacity. This shift reveals a key system dynamic: when prevented from oversizing generation, the model compensates with massive storage deployment. However, with lithium-ion batteries limited to 6 hours of storage duration, this solution may struggle to bridge extended periods of low renewable generation, suggesting the need for longer-duration storage technologies.

Baseline with renewable generation limits and flow battery scenario#

In this scenario, in addition to the upper limit imposed on renewable energy generation we also assume also flow batteries are available for energy storage.

n = create_network(ts, costs, add_flow_battery=True, limit_renewable_energy=True)

n.optimize(
    solver_name="highs",
    solver_options=SOLVER_OPTIONS,
)

Hide code cell output

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
/home/thesethtruth/.cache/pypoetry/virtualenvs/esi-course-6syH-QiP-py3.12/lib/python3.12/site-packages/linopy/common.py:154: UserWarning: coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.
  warn(
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options:
 - primal_feasibility_tolerance: 1e-05
 - dual_feasibility_tolerance: 1e-05
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 16/16 [00:00<00:00, 27.36it/s]
Writing continuous variables.: 100%|██████████| 6/6 [00:00<00:00, 62.17it/s]
INFO:linopy.io: Writing time: 0.7s
Running HiGHS 1.8.1 (git hash: 4a7f24a): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-03, 1e+02]
  Cost   [5e-01, 3e+03]
  Bound  [0e+00, 0e+00]
  RHS    [2e+04, 3e+05]
Presolving model
91953 rows, 65677 cols, 249579 nonzeros  0s
91953 rows, 65677 cols, 249579 nonzeros  0s
Presolve : Reductions: rows 91953(-74493); columns 65677(-4407); elements 249579(-78900)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 8760(4.83545e+08) 0s
      12160     1.8623553094e+08 Pr: 14874(3.04767e+10); Du: 0(0.00164249) 6s
      17198     1.8718267060e+08 Pr: 17058(7.99661e+10); Du: 0(0.00204629) 11s
      21942     1.8843983537e+08 Pr: 19082(1.40572e+11); Du: 0(0.00272412) 17s
      25860     1.8951551343e+08 Pr: 20686(1.90755e+11); Du: 0(0.00344047) 23s
      28739     3.1744926814e+08 Pr: 21302(2.66921e+11); Du: 0(0.00379979) 31s
      31281     3.6779305938e+08 Pr: 21078(7.42084e+10); Du: 0(0.00353308) 36s
      34243     5.8487688436e+08 Pr: 37451(9.3575e+11); Du: 0(0.00281197) 42s
      38652     6.3551931227e+08 Pr: 17331(1.10048e+11); Du: 0(0.00264149) 47s
      45183     9.2581428389e+08 Pr: 36514(2.68435e+10); Du: 0(0.00390185) 54s
      54384     1.0985648279e+09 Pr: 14907(2.37384e+09); Du: 0(0.00450506) 60s
      64289     1.1577672564e+09 Pr: 3050(1.17193e+08); Du: 0(0.00353224) 66s
      76312     1.1844667613e+09 Pr: 0(0); Du: 0(7.02838e-12) 70s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-84l_bh4r
Model status        : Optimal
Simplex   iterations: 76312
Objective value     :  1.1844667613e+09
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 70084 primals, 166446 duals
Objective: 1.18e+09
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
Relative P-D gap    :  7.4476445564e-15
HiGHS run time      :         69.70
Writing the solution to /tmp/linopy-solve-afoddizw.sol
('ok', 'optimal')

Plot the energy system configuration#

fig, ax = plt.subplots(1, 1, figsize=(10, 3))
plt.subplots_adjust(hspace=0.4)

capacity_data = pd.concat(
    [
        n.generators.p_nom_opt,
        n.storage_units.p_nom_opt,
    ],
    axis=1,
).div(1e3)

capacity_data.columns = ["generation", "storage"]

capacity_data.plot(kind="barh", ax=ax)
ax.set_title("Installed Capacities: Baseline with LDES and limited renewable energy")
ax.set_ylabel("Capacity (GW)")
Text(0, 0.5, 'Capacity (GW)')
_images/e73de76336a07e21e709d73b80bd8264aa5603c031a69d922ad12b05de9276cd.png

Impact of adding long-duration energy storage#

After adding flow batteries (100-hour storage duration), the model selects:

  • 300 GW solar (at limit)

  • 150 GW wind (at limit)

  • 80 GW battery electric storage

  • 40 GW flow battery storage

The introduction of long-duration energy storage (LDES) reduces the need for short-duration batteries by about 80%. The model opts for a mix of both storage types, using flow batteries to manage multi-day supply gaps while keeping some battery capacity for daily cycling. This more diversified storage portfolio better matches the temporal patterns of renewable generation.

Sensitivity analysis#

Define scenarios and run the optimisation problems#

# Create and solve networks with different cost scenarios
scenarios = {
    "baseline": (1.0, 1.0),
    "expensive_generation": (2.0, 1.0),
    "expensive_storage": (1.0, 2.0),
}

results = {}
for name, (gen_mult, store_mult) in scenarios.items():
    n = create_network(
        ts,
        costs,
        gen_mult,
        store_mult,
        add_flow_battery=True,
        limit_renewable_energy=True,
    )
    n.optimize(solver_name="highs", solver_options=SOLVER_OPTIONS)
    results[name] = n

Hide code cell output

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
/home/thesethtruth/.cache/pypoetry/virtualenvs/esi-course-6syH-QiP-py3.12/lib/python3.12/site-packages/linopy/common.py:154: UserWarning: coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.
  warn(
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options:
 - primal_feasibility_tolerance: 1e-05
 - dual_feasibility_tolerance: 1e-05
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 16/16 [00:00<00:00, 27.89it/s]
Writing continuous variables.: 100%|██████████| 6/6 [00:00<00:00, 59.46it/s]
INFO:linopy.io: Writing time: 0.7s
Running HiGHS 1.8.1 (git hash: 4a7f24a): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-03, 1e+02]
  Cost   [5e-01, 3e+03]
  Bound  [0e+00, 0e+00]
  RHS    [2e+04, 3e+05]
Presolving model
91953 rows, 65677 cols, 249579 nonzeros  0s
91953 rows, 65677 cols, 249579 nonzeros  0s
Presolve : Reductions: rows 91953(-74493); columns 65677(-4407); elements 249579(-78900)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 8760(4.83545e+08) 0s
      12160     1.8623553094e+08 Pr: 14874(3.04767e+10); Du: 0(0.00164249) 6s
      18075     1.8738298523e+08 Pr: 17422(9.16745e+10); Du: 0(0.00227149) 11s
      23161     1.8881790897e+08 Pr: 19606(1.58739e+11); Du: 0(0.00280235) 17s
      27088     1.8986007861e+08 Pr: 21214(2.07089e+11); Du: 0(0.00374331) 23s
      28739     3.1744926814e+08 Pr: 21302(2.66921e+11); Du: 0(0.00379979) 29s
      31281     3.6779305938e+08 Pr: 21078(7.42084e+10); Du: 0(0.00353308) 35s
      34243     5.8487688436e+08 Pr: 37451(9.3575e+11); Du: 0(0.00281197) 40s
      38652     6.3551931227e+08 Pr: 17331(1.10048e+11); Du: 0(0.00264149) 46s
      44248     8.2510192110e+08 Pr: 17787(1.64779e+10); Du: 0(0.00397365) 52s
      49751     1.0131698839e+09 Pr: 32853(5.91491e+09); Du: 0(0.00406809) 57s
      58469     1.1310300563e+09 Pr: 5691(4.20632e+08); Du: 0(0.00401105) 65s
      72840     1.1814921064e+09 Pr: 1640(1.93166e+07); Du: 0(0.00526611) 71s
      76312     1.1844667613e+09 Pr: 0(0); Du: 0(7.02838e-12) 72s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-2jquwdks
Model status        : Optimal
Simplex   iterations: 76312
Objective value     :  1.1844667613e+09
Relative P-D gap    :  7.4476445564e-15
HiGHS run time      :         71.85
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 70084 primals, 166446 duals
Objective: 1.18e+09
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
Writing the solution to /tmp/linopy-solve-igi72umu.sol
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
/home/thesethtruth/.cache/pypoetry/virtualenvs/esi-course-6syH-QiP-py3.12/lib/python3.12/site-packages/linopy/common.py:154: UserWarning: coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.
  warn(
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options:
 - primal_feasibility_tolerance: 1e-05
 - dual_feasibility_tolerance: 1e-05
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 16/16 [00:00<00:00, 28.02it/s]
Writing continuous variables.: 100%|██████████| 6/6 [00:00<00:00, 66.72it/s]
INFO:linopy.io: Writing time: 0.69s
Running HiGHS 1.8.1 (git hash: 4a7f24a): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-03, 1e+02]
  Cost   [5e-01, 3e+03]
  Bound  [0e+00, 0e+00]
  RHS    [2e+04, 3e+05]
Presolving model
91953 rows, 65677 cols, 249579 nonzeros  0s
91953 rows, 65677 cols, 249579 nonzeros  0s
Presolve : Reductions: rows 91953(-74493); columns 65677(-4407); elements 249579(-78900)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 8760(4.83545e+08) 0s
      12154     1.8629435974e+08 Pr: 14919(2.6928e+10); Du: 0(0.00161685) 6s
      17486     1.8747845492e+08 Pr: 17240(8.07018e+10); Du: 0(0.00220029) 11s
      22167     1.8899507311e+08 Pr: 19197(1.41553e+11); Du: 0(0.00276781) 16s
      26578     1.9089992278e+08 Pr: 21124(2.26527e+11); Du: 0(0.00285984) 26s
      29890     3.4561171789e+08 Pr: 21653(1.36529e+11); Du: 0(0.00317713) 32s
      33581     3.7300537692e+08 Pr: 21521(1.98926e+11); Du: 0(0.0026611) 37s
      35501     6.8653366632e+08 Pr: 17395(2.70996e+10); Du: 0(0.00262549) 43s
      38459     9.2146333615e+08 Pr: 34665(2.76269e+10); Du: 0(0.00255441) 49s
      43425     1.0379168700e+09 Pr: 34471(4.58596e+10); Du: 0(0.00336107) 56s
      45541     1.0725731154e+09 Pr: 34647(4.80194e+10); Du: 0(0.00346386) 62s
      51089     1.2588116694e+09 Pr: 12861(2.38788e+09); Du: 0(0.00310731) 68s
      56524     1.3764401555e+09 Pr: 7566(1.2695e+09); Du: 0(0.00308137) 73s
      62344     1.4633140902e+09 Pr: 2423(4.93041e+07); Du: 0(0.00439701) 79s
      64597     1.4694272694e+09 Pr: 0(0); Du: 0(3.39427e-11) 81s
      64597     1.4694272694e+09 Pr: 0(0); Du: 0(3.39427e-11) 81s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-1nx5v352
Model status        : Optimal
Simplex   iterations: 64597
Objective value     :  1.4694272694e+09
Relative P-D gap    :  5.1920872098e-15
HiGHS run time      :         81.12
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 70084 primals, 166446 duals
Objective: 1.47e+09
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.
Writing the solution to /tmp/linopy-solve-i6qfaceh.sol
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity'], dtype='object', name='Bus')
/home/thesethtruth/.cache/pypoetry/virtualenvs/esi-course-6syH-QiP-py3.12/lib/python3.12/site-packages/linopy/common.py:154: UserWarning: coords for dimension(s) ['Generator'] is not aligned with the pandas object. Previously, the indexes of the pandas were ignored and overwritten in these cases. Now, the pandas object's coordinates are taken considered for alignment.
  warn(
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options:
 - primal_feasibility_tolerance: 1e-05
 - dual_feasibility_tolerance: 1e-05
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 16/16 [00:00<00:00, 26.07it/s]
Writing continuous variables.: 100%|██████████| 6/6 [00:00<00:00, 58.14it/s]
INFO:linopy.io: Writing time: 0.74s
Running HiGHS 1.8.1 (git hash: 4a7f24a): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-03, 1e+02]
  Cost   [5e-01, 4e+03]
  Bound  [0e+00, 0e+00]
  RHS    [2e+04, 3e+05]
Presolving model
91953 rows, 65677 cols, 249579 nonzeros  0s
91953 rows, 65677 cols, 249579 nonzeros  0s
Presolve : Reductions: rows 91953(-74493); columns 65677(-4407); elements 249579(-78900)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 8760(4.83545e+08) 0s
      12703     2.0489470433e+08 Pr: 13108(6.74711e+10); Du: 0(0.000110946) 6s
      16219     2.0639473589e+08 Pr: 14226(1.25799e+11); Du: 0(0.000113852) 12s
      20599     2.0876460028e+08 Pr: 15557(2.17646e+11); Du: 0(6.5762e-05) 17s
      24885     2.1180214621e+08 Pr: 16805(3.37596e+11); Du: 0(2.87633e-05) 22s
      28780     3.6006370339e+08 Pr: 22921(8.00132e+09); Du: 0(0.000254273) 28s
      32743     7.0655272102e+08 Pr: 35444(5.189e+10); Du: 0(0.00108661) 35s
      37263     8.3214993063e+08 Pr: 35309(2.11233e+10); Du: 0(0.0015753) 41s
      39920     8.5589755519e+08 Pr: 20430(1.62566e+10); Du: 0(0.00211412) 50s
      42957     9.1996282976e+08 Pr: 34670(2.82353e+10); Du: 0(0.00245885) 56s
      46898     9.7758558247e+08 Pr: 17530(7.7805e+09); Du: 0(0.0029493) 61s
      51334     1.0595044871e+09 Pr: 14641(9.76205e+09); Du: 0(0.00327851) 66s
      56623     1.1927766045e+09 Pr: 33264(1.46378e+10); Du: 0(0.00323071) 72s
      62079     1.2694809687e+09 Pr: 14228(4.87321e+09); Du: 0(0.00450072) 78s
      64547     1.3397520990e+09 Pr: 12356(2.37925e+09); Du: 0(0.00452715) 84s
      69694     1.4799065430e+09 Pr: 36689(8.19745e+09); Du: 0(0.00480656) 90s
      73479     1.5311009901e+09 Pr: 5398(3.80298e+08); Du: 0(0.00532698) 95s
      78428     1.5646680992e+09 Pr: 4826(1.59483e+09); Du: 0(0.00574353) 102s
      79839     1.5662568840e+09 Pr: 0(0); Du: 0(1.7585e-10) 103s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-1w46vs78
Model status        : Optimal
Simplex   iterations: 79839
Objective value     :  1.5662568840e+09
Relative P-D gap    :  8.9810913591e-15
HiGHS run time      :        103.64
Writing the solution to /tmp/linopy-solve-e108d6c5.sol
INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 70084 primals, 166446 duals
Objective: 1.57e+09
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.

Results visualisation#

Let’s visualise the results so we can compare the results from our different scenarios.

fig, axes = plt.subplots(1, 2, figsize=(12, 5))
plt.subplots_adjust(hspace=0.4)

# Capacity comparison
capacity_data = pd.DataFrame()
for scenario, n in results.items():
    caps = pd.concat([n.generators.p_nom_opt, n.storage_units.p_nom_opt]).rename(
        scenario
    )
    capacity_data = pd.concat([capacity_data, caps], axis=1)

# convert to GW
capacity_data = capacity_data.div(1e3)
# make scenario names pretty
capacity_data.columns = [name.replace("_", " ").capitalize() for name in capacity_data.columns]

capacity_data.plot(kind="barh", ax=axes[0])
axes[0].set_title("Installed Capacities by scenario")
axes[0].set_xlabel("Capacity (GW)")

# Add data call outs
for container in axes[0].containers:
    axes[0].bar_label(container, fmt='%.0f', label_type='edge', fontsize=8, padding=3)

# Cost comparison
cost_data = pd.DataFrame(
    {name: n.objective for name, n in results.items()}, index=["Relative System Cost"]
).T.div(1e9)

# Normalize to baseline
cost_data = cost_data.div(cost_data.loc["baseline"]).mul(100)
# make index pretty
cost_data.index = [name.replace("_", " ").capitalize() for name in cost_data.index]

cost_data.plot(kind="barh", ax=axes[1], legend=False)
axes[1].set_title("Relative system cost by scenario")
axes[1].set_xlabel("Relative cost (%)")

for container in axes[1].containers:
    axes[1].bar_label(container, fmt="%.0f%%", label_type="center", fontsize=12, color="white", fontweight="bold")

plt.tight_layout()
_images/9c8e411f1a32c4cf0ec5fab781a44b1986335e4efc53ef18500ccc32c16c4502.png

Sensitivity analysis results#

Comparing three cost scenarios reveals clear system adaptation patterns:

  • In the baseline scenario (renewable and storage costs at reference values), the model deploys 150 GW wind, 300 GW solar, with a mixed storage portfolio of 51 GW batteries and 44 GW flow batteries.

  • When generation costs double, wind capacity slightly decreases to 143 GW while solar remains at its 300 GW limit. Storage capacities remain similar (50 GW batteries, 50 GW flow), suggesting storage deployment is relatively insensitive to generation costs. Total system costs increase by 30%.

  • With doubled storage costs, the model shifts notably. While solar remains capped and wind drops slightly to 140 GW, battery storage is completely phased out in favor of increased flow battery capacity (83 GW). This suggests flow batteries’ longer duration becomes more valuable when storage is expensive. The scenario results in the highest cost increase at 34% above baseline. Interestingly, the deployed wind capacity is more sensitive towards the costs of storage than towards the cost of generation.

Detailed analysis of system dynamics#

We might be interested in specific system dynamics of one of the scenarios during a specific week. Remember the first chapter of this tutorial? We found that February has a long period of high demand due to heating. Let’s see how each of the energy system configurations fares during this period.

def plot_timeseries(network, scenario_name, week_start=0):
    week_slice = slice(week_start * 168, (week_start + 1) * 168)

    fig, ax = plt.subplots(figsize=(15, 6))

    # slice the generation data
    generation = network.generators_t.p_max_pu[week_slice] * network.generators.p_nom_opt
    storage = network.storage_units_t.p[week_slice]
    load = network.loads_t.p[week_slice]

    # extract the negative and positive parts of the storage
    storage_negative = storage.copy()
    storage_negative[storage > 0] = 0
    storage_negative.columns = [
        f"{col} (discharging)" for col in storage_negative.columns
    ]

    storage_positive = storage.copy()
    storage_positive[storage < 0] = 0
    storage_positive.columns = [f"{col} (charging)" for col in storage_positive.columns]

    # combine the parts to positive and negative timeseries (on the bus)
    positives = pd.concat([generation, storage_positive], axis=1).div(1e3)
    negatives = pd.concat([-load, storage_negative], axis=1).div(1e3)

    curtailment = -(positives.sum(axis=1) + negatives.sum(axis=1))
    curtailment[curtailment > 0] = 0
    curtailment.name = "curtailment"

    negatives = pd.concat([negatives, curtailment], axis=1)

    # Define colors for each technology
    colors = {
        "wind": "#4B90B9",
        "solar": "#E6B829",
        "demand": "#708090",
        "battery storage (charging)": "#9CAF88",
        "battery storage (discharging)": "#4A5D3F",
        "flow battery (charging)": "#9B8EA9",
        "flow battery (discharging)": "#563D67",
        "curtailment": "#B07979",
    }
    # Apply colors to the plots
    positives.plot(
        kind="area",
        ax=ax,
        stacked=True,
        color=[colors.get(col, "#333333") for col in positives.columns],
    )
    negatives.plot(
        kind="area",
        ax=ax,
        stacked=True,
        color=[colors.get(col, "#333333") for col in negatives.columns],
    )

    ax.set_ylim(-250, 250)
    ax.set_title(f"Power System Operation - Week {week_start+1} - {scenario_name}")
    ax.set_ylabel("Power (GW)")
    ax.set_xlabel("")
    plt.legend(bbox_to_anchor=(1.05, 1), loc="upper left", title="Technology", frameon=False)

    return fig


fig = plot_timeseries(
    results["baseline"], "Baseline", week_start=4
)
fig = plot_timeseries(
    results["expensive_generation"], "Expensive Generation", week_start=4
)
fig = plot_timeseries(
    results["expensive_storage"], "Expensive Storage", week_start=4
)
/home/thesethtruth/.cache/pypoetry/virtualenvs/esi-course-6syH-QiP-py3.12/lib/python3.12/site-packages/pandas/plotting/_matplotlib/core.py:1800: UserWarning: Attempting to set identical low and high ylims makes transformation singular; automatically expanding.
  ax.set_ylim(None, 0)
/home/thesethtruth/.cache/pypoetry/virtualenvs/esi-course-6syH-QiP-py3.12/lib/python3.12/site-packages/pandas/plotting/_matplotlib/core.py:1800: UserWarning: Attempting to set identical low and high ylims makes transformation singular; automatically expanding.
  ax.set_ylim(None, 0)
/home/thesethtruth/.cache/pypoetry/virtualenvs/esi-course-6syH-QiP-py3.12/lib/python3.12/site-packages/pandas/plotting/_matplotlib/core.py:1800: UserWarning: Attempting to set identical low and high ylims makes transformation singular; automatically expanding.
  ax.set_ylim(None, 0)
_images/d33b56d0d756e49b4befb9e4dfd015fdc429de38633a294b862b6e54bab2f9df.png _images/e89af14bad7b2e1bf60d8cafe8750f266576f16b08b83463fb23099c9cbec82c.png _images/65a1412a5f88574fab0c43ccea5178824a0b5067a27b207160f4e6e788cc6e9e.png

Conclusion#

  1. The addition of long-duration storage (flow battery) allows for better management of longer periods of low renewable generation.

  2. Different cost scenarios lead to different optimal system configurations:

    • Higher generation costs result in more storage capacity to make better use of available renewable energy

    • Higher storage costs push the system toward oversized generation capacity to reduce storage needs

  3. The baseline scenario represents a balanced approach between generation and storage.

Next steps#

To extend this analysis, you could:

  1. Add more technology options (e.g., different storage or dispatchables technologies)

  2. Analyze seasonal patterns using the full year of data

  3. Include more detailed cost components (fixed O&M, variable O&M, etc.)